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ABSTRACT 


L  cxY\\  'cla. 

~4'his  thesis  presents  an  analytical^sfrid^  of  laminar  convective  heat  transfer  along  vertical 


slender  cylinders,  Chapter  2  covers  the  case  of  variable  surface  heating  conditions  ol  the  form 

cf  ■  V  uj  '  c- y  i  .  i  i  4  V  v  >■  ,_~T 

tC(x)  =(TOn)for  natural  convection.  Results  are  presented  for  Pr  =  0.1,  0.7,  7,  and  100  and  for  n 


=  -0.5,  -0.25,  0,  0.25,  and  0.5.  It  is  found  that  as  the  curvature  parameter  A)  increases  the  local 
heat  transfer  coefficient  increases  for  the  range  of  n  investigated.  The  local  wall  shear  stress 


decreases  for  a  given  location,*  with  increasing  curvature.  It  increases  with  increasing  x  for  a 
constant  radius  cylinder  with  rt^O,  but  may  decrease  for  n  <  0.  It  is  also  found  that  a  higher 


value  of  n  or  Prandtl  number  resulted  in  a  lower  wall  velocity  gradient,  and  a  higher  value  for 


Nu,GrV  s-  llowever,  as  the  curvature  is  increased,  the  effect  of  n  and  Pr  on  these  values 
diminishes. 


Chapter  3  covers  the  entire  range  of  buoyancy-assisted,  mixed  convection  along  a  vertical 
slender  cylinder  for  the  case  of  variable  wall  temperature  distribution.^  of  the  form 
I  w(x)  -  I  „  =  uxn.  Vhe  buoyancy  parameter  varied  from  zero  for  pure  free  convection  to 


one  for  pure  forced  convection. ^  Results  are  presented  for  Pr  =  0.1,  0.7,  7,  and  100  and  for  n  - 

ej .  .  <£Ep 

-0.4,  -0,2,  0,  0.2  and  0.5.  Tt  is  found  that  for  given  values  of  Pr,  n  and  f,,  as  the  curvature 
*■ — L  _ci  l,  do- 

parameter^  increases  both  the  local  wall  shear  stress  and  the  local  heat  transfer  coefficient 


increase.  higher  values  of  n  and  Pr  are  found  to  lower  the  local  wall  shear  stress  and  to 

raise  the  local  Nusselt  number,  and  the  effect  of  n  and  Pr  is  found  to  lessen  as  the  curvature 
parameter  is  increased^  Average  Nusselt  numbers  were  calculated  and  correlation  equations  for 
the  local  a/id  average  Nusselt  numbers  are  presented  for  for  both  problems  investigated 
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radial  coordinate 
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thermal  diffusivity 
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curvature  parameter  for  natural  convection  with  variable  surface  temperature, 
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kinematic  viscosity 

mixed  convection  parameter  for  variable  surface  temperature,  (1  t  UJ ‘)  1 
mixed  convection  parameter  for  variable  surface  heat  flux,  (1  f  ')  1 
buoyancy  parameter  for  variable  surface  temperature,  (ir,/Re* 
buoyancy  parameter  for  variable  surface  heat  flux,  (Jr*, /Re* 2 


>\v 
‘‘  V 

L  O  - 

V  v 

vv 

V  .v  . 


1 


vV 

•VV 

-\\v> 


i 


rm 


/■.S.s 

y.>y. 

rv  vs 


•  •  «  M  A 


V  V  V  ' 

•v-y-/. 

.’'Vv'vj 


:;y:^ 


a o  v: < ’.-.y. 


v,v 

y 

A  ^ .  »*. 
•/V  V  V 


INTRODUCTION 


This  thesis  studies  convective  heat  transfer  along  slender  vertical  cylinders.  Chapter  2 
presents  an  analysis  of  natural  convection  with  a  power  law  variation  in  surface  heat  flux  of  the 
form  q*(x)  =  ax".  In  Chapter  3  mixed  convection  covering  the  full  range  from  pure  free 
convection  to  pure  forced  convection  for  a  power  law  variation  in  surface  temperature  of  the 
form  Tw(x)  -  T *  =  ax"  is  studied. 

In  formulating  the  problem  of  free  convection  along  slender  vertical  cylinders,  the 
governing  equations  are  written  in  their  cylindrical  coordinate  form  and  then  transformed  into  a 
dimensionless  form.  The  boundary  layers  are  nonsimilar  due  to  the  transverse  surface  curvature. 
As  the  radius  rn  approaches  infinity,  the  equations  reduce  to  those  for  a  vertical  plate,  and  the 
boundary  layers  become  similar.  The  vertical  plate  problem  has  been  studied  extensively.  The 
isothermal  vertical  flat  plate  was  first  analyzed  by  Ostrach  [1].  Sparrow  and  Gregg  [2]  later 
studied  the  case  of  uniform  surface  heat  flux.  Subsequently  many  analytical  studies  using 
various  solution  methods  were  used  to  solve  variations  of  the  problem  (see,  e.g.  Kefs.  [3  -  7]). 
These  studies  included  cases  of  variable  surface  temperature  distributions.  The  case  of  the 
power  law  temperature  variation  or  a  power  law  variation  in  surface  heat  flux  lends  itself  to 
similarity  solutions  for  the  vertical  flat  plate.  Results  for  these  can  be  found  in  Spturow  and 
Gregg  [5]  for  the  power  law  variation  in  surface  temperature  and  in  Chen  et  al.  [7]  for  both 
cases  that  include  the  local  nonsimilarity  solution  for  inclined  plates. 

The  case  of  pure  forced  convection  past  a  flat  plate  is  the  well  known  Blasius  solution  for 
the  momentum  equation  and  the  Pohlhausen  solution  for  the  case  of  uniform  wall  temperature 
(LAVT).  The  case  of  uniform  surface  heat  flux  (CUT)  was  solved  numerically  for  the  entire 
range  of  Prandtl  numbers  (Pr)  by  Churchill  and  O/.oe  [8].  The  solutions  for  both  pure  free 
convection  and  pure  forced  convection  have  been  accurately  correlated  and  should  provide  a 
basis  for  comparison  for  the  limit  of  the  vertical  cylinder  problem  as  r„  -*  ^o. 

The  case  of  pure  forced  convection  past  a  vertical  cylinder  has  been  studied  by  various 
authors  [9—  14],  Seban  and  Bond  [9]  conducted  an  early  study  considering  the  l  AVI'  case  for 
Pr  =  0.715  and  curvatures  (A,.  =  2  —Re,1'2)  up  to  3.0  using  a  power  series  method  of  solution. 

rn 
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Chen  and  Mucoglu  [13]  first  presented  results  lor  the  case  of  UWT  as  part  of  the  solution  for 
mixed  convection  and  then  fallowed  that  with  the  case  of  (Jill  [14],  Their  results  covered 
curvatures  up  to  AF  =  8  for  Pr  =  0.7  and  7,  and  they  used  the  local  nonsimilanty  solution.  Bui 
and  Cebeci  [12]  analyzed  the  UWT  case  as  part  of  their  mixed  convection  study  for  Pr  ~  0.1, 
1.0,  and  10  lor  curvatures  up  to  AF  =  10  using  a  central  difference  finite-difference  method  of 
solution.  Ice  et  a!.  [10,  11]  studied  both  the  UWT  and  Ulli  cases  as  part  of  their  mixed 
convection  analysis  for  Pr  =  0.1,  0.7,  7,  and  100  for  curvature  parameter  A,,  up  to  50  using  a 
weighted  finite-difference  scheme.  Tee  et  al.  qualified  the  accuracy  of  their  results  for  A,  >  10 
due  to  the  use  of  an  insufficiently  small  step  size  in  the  radial  direction. 

More  extensively  researched  has  been  the  case  of  natural  convection  along  a  cylinder 
[10,11,15-22].  I-lcnbaas  [15]  used  Uangmuir's  stagnant  film  model  to  evaluate  the  heat 
transfer  coefficients  for  vertical  cylinders  under  the  UWT  condition.  Sparrow  and  Gregg  [16] 
refined  this  method  and  reworked  this  problem  using  a  power  series  solution  and  obtained 
results  for  Pr  =  0.72  and  1  covering  curvature  parameter,  AN  =  2— - tir, 1  \  up  to  about  1.5. 

*o 

Kiitkcn  [IS]  also  used  a  power  series  expansion  and  investigated  power  law  variation  in  wall 
temperature  and  the  I  MP  condition  for  0.7  <  Pr  <  10.  I  ujii  and  Uehara  [20]  also  applied  a 
power  series  solution  to  cylinders  with  variable  wall  temperature  for  AN  up  to  2,73  and  for 
0.7  <  Pr  <  100.  To  reduce  the  truncation  errors  of  the  power  series  solution  and  the  local 
similarity  solution,  Minkowycz  and  Sparrow  [21]  employed  the  local  nonsimilarity  solution 
method  and  obtained  results  for  Pr  =  0.733  and  0  <  AN  <  10.  I,ee  et  al.  [10,  11]  solved  the 
UWT  and  the  l,’ IIP  cases  of  natural  convection  along  a  vertical  cylinder  for  Pr  =  0.1,  0.7,  7, 
and  100  as  part  of  their  mixed  convection  study  employing  a  finite-difference  method  designed 
to  overcome  the  difficulties  associated  with  increasing  curvature.  As  a  result,  I.ee  et  al  carried 
their  results  to  AN  =  50.  More  recently,  I.ee  et  al.  [17]  used  a  modified  method  [23]  to  solve 
the  variable  wall  temperature  problem  for  natural  convection  along  a  vertical  cylinder.  They 
found  that  the  results  of  [10]  were  incorrect  for  AN  >  10  due  to  an  improper  step  size  used  in 
the  radial  direction  . 
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The  significant  finding  for  pure  forced  or  pure  free  convection  along  a  vertical  cylinder  is 
that  the  nonsimilarity  due  to  the  surface  curvature  increases  the  heat  transfer  coefficient.  The 
velocity  gradient  at  the  wall  increases  with  increasing  surface  curvature  for  pure  free  convection 
under  a  power  law  variation  in  surface  temperature  or  for  pure  forced  convection,  but  decreases 
for  the  UHF  case  for  free  convection.  Minkowycz  and  Sparrow  [21]  showed  that  the  error 
caused  by  the  power  series  solution  increases  with  increasing  curvature.  Ixe  et  al.  [17]  showed 
that  the  results  provided  by  Fujii  and  L'ehara  [20]  are  valid  only  up  to  As  of  about  1.5.  Fee  et 
al.  [10]  demonstrated  that  the  central  difference  scheme  used  by  Bui  and  Cebeci  [12]  produces 
large  errors  as  the  curvature  increases,  but  the  results  in  [10,  1 1]  are  not  accurate  for  A  >  10 
because  they  used  too  large  a  step  size  in  the  radial  direction. 

There  have  been  relatively  few  studies  on  mixed  convection  along  a  vertical  cylinder.  The 

studies  by  Bui  and  Cebeci  [12]  and  Chen  and  Mucoglu  [13]  are  for  the  UWT  case.  They 

investigated  the  forced  convection  dominated  regime  in  terms  of  the  buoyancy  force  parameter 

£ix  =  Grx/Re].  Mucoglu  and  Chen  also  presented  the  UI1F  case  with  the  modified  buoyancy 

parameter  Ci*x  —  Gr*x/Re^/2.  Lee  et  al.  [10]  treated  the  UWT  case  for  the  entire  regime  using  a 

new  mixed  convection  parameter,  y  =  ( I  4-  £2J/4)" which  varies  from  zero  for  pure  free 

convection  to  one  for  pure  forced  convection.  In  Lee  et  al.  [11],  the  UHF  case  is  presented  for 

the  entire  regime  using  y*  =  (1  +  12+i/s)1.  The  main  difficulty  addressed  in  these  two  studies  has 

been  the  effect  of  curvature  on  the  numerical  solution.  T  he  results  of  Chen  and  Mucoglu 

[13,  14]  for  the  forced-convection  dominated  solutions  agree  well  with  those  of  [10,  11]  for 

curvature  parameters  [A  =  2-j^-(ReJ'2  +  Gr[/4)  1  or  A*  =  2-£-(ReJ/2  +  Gr*[/S)  ']  from  0  to  X  for 

*0  *0 

Pr  =  0.7  and  7.  T  he  results  of  Bui  and  Cebeci  for  large  curvature  parameters  were  found  to  be 
inaccurate  by  Lee  et  al.  [10]  due  to  their  use  of  a  central  difference  finite-difference  method  of 
solution.  1-ee  et  al.  [10,  1 1]  used  a  weighted  finite-difference  scheme  to  correct  this  but  failed  to 
use  the  proper  step  size  for  large  curvatures. 

For  comparison  purposes  the  mixed  convection  solution  for  the  vertical  flat  plate  is 
needed.  There  are  numerous  studies  of  this  flow  geometry  (sec,  e.g.  Refs.  [24—  20]).  T  hese 
numerical  results  have  been  correlated  by  Chen  et  al.  [24]  and  Churchill  [30]  using  the  simple 
form  Num  =  \uFm  +  NuNm  with  m  =  3,  in  which  NuF  is  the  Nussclt  number  for  pure  forced 


convection  and  Nun  is  the  Nusselt  number  for  pure  free  convection.  Ibe  correlation  was  also 
verified  experimentally  for  the  case  of  air  [24]. 

To  date  no  study  has  been  done  for  the  case  of  variable  surface  heat  flux  along  a  slender 
vertical  cylinder  for  natural  convection.  Ibis  will  be  undertaken  in  Chapter  2.  Similarly,  the 
combined  effects  of  buoyancy  force  and  curvature  in  mixed  convection  with  variable  surface 
temperature  along  a  slender  vertical  cylinder  has  not  been  solved,  furthermore,  the  results  of 
Lee  et  al.  [10]  and  Bui  and  Cebeci  [12]  for  the  UWT  case  are  questionable  for  high  values  of 
the  curvature  parameter  and  need  to  be  reassessed.  'Ibis  is  the  topic  of  Chapter  3. 
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NATURAL  CONVECTION  ALONG  A  VERTICAL  CYLINDER 
WITH  VARIABLE  SURFACE  HEAT  FLUX 


J.  J.  Heckel,  T.  S.  Chen,  and  B.  F.  Armaly 
Department  of  Mechanical  and  Aerospace  Engineering 
University  of  Missouri-Rolla,  Rolla,  MO  65401 


abstract 


Natural  convection  in  laminar  boundary  layer  flow  along  slender  vertical  cylinders  is 
analyzed  for  the  situation  in  which  the  surface  heat  flux  qw(x)  varies  arbitrarily  with  the  axial 
coordinate  x.  The  governing  boundary  layer  equations  along  with  the  boundary  conditions  are 
first  cast  into  a  dimensionless  form  by  a  nonsimilar  transformation  and  the  resulting  system  of 
equations  is  then  solved  by  a  weighted  finite-difference  method  of  solution  in  conjunction  with 
the  cubic  spline  interpolation.  Sample  calculations  were  performed  for  the  case  of  power  law 
variation  in  surface  heat  flux,  q„(x)  =  ax",  for  fluids  with  Prandtl  numbers  of  0.1,  0.7,  7,  and  100 
over  a  wide  range  of  values  for  the  curvature  parameter  0  <,  A+N  <,  50  or  0  <,  <J  <  5.  Results  for 
the  local  and  average  Nusselt  numbers  as  well  as  the  local  wall  shear  stress  were  obtained.  It  is 
found  that  the  value  of  the  radial  coordinate  rj^  must  be  increased,  as  the  curvature  parameter  is 
increased,  to  obtain  valid  results.  Also  the  local  Nusselt  number  is  found  to  increase  with 
increasing  curvature,  Prandtl  number,  and  the  exponent  n.  For  any  given  location  of  x,  the 
local  wall  shear  stress  was  found  to  decrease  with  increasing  curvature,  increasing  Prandtl 
number,  and  increasing  n.  Correlation  equations  for  the  local  and  average  Nusselt  numbers  are 
presented. 
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GrY,  Gr*0 
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Nu, 
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nomi:nci.ahjri: 

reduced  stream  function,  (^/5vr„)(Gr//5)  1/5 

<?f /d»|  or  gravitational  acceleration 

local  Grashof  number,  g/?[T„(x)  -  l’00]xJ/v2 

modified  local  Grashof  number,  g/?qw(x)x4/A:v2 

modified  Grashof  numbers  defined,  respectively,  as 
g/?qw(I.)I  4/^v2  and  gflqw(Toyjkv2 

local  heat  transfer  coefficient,  qw/(Tw  — 

average  heat  transfer  coefficient 

thermal  conductivity  of  the  fluid 

an  arbitrary  length  of  cylinder 

exponent  in  the  power  law  variation  of  heat  flux 

local  Nusselt  number,  hxjk 

average  Nusselt  number,  hi. Ik 

Prandtl  number,  v/a 

radial  coordinate 

radius  of  the  cylinder 

fluid  temperature 

axial  velocity  component 

radial  velocity  component 

axial  coordinate 


(I 

Y 
* 1 
0 
X 

A,^ 

AY 


thermal  diffusivity 

volumetric  coefficient  of  thermal  expansion 
variable  heat  flux  parameter,  (x/qw)(dqw/dx) 
pseudo-similarity  variable,  [(r2  -  r2)/2r0x](Gr//5)l/s 
dimensionless  temperature,  (T  —  I ’„)(Gr,f/5)1/5/[qw(x)x/<:] 
curvature  parameter,  (2x/r„)(Gr,*/5)  1/5 

curvature  parameter  for  variable  surface  temperature,  (2x/r0)(  irx  1 4 
modified  curvature  parameter,  (2x/r„)Gr\ 1,5 
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dynamic  viscosity 
kinematic  viscosity 

dimensionless  axial  coordinate,  (A*N/2),/2 

fluid  density 

local  wall  shear  stress 

normalized  temperature  profile,  <£(£,  >;)  =  0({,  >?)/0(£,  0) 
stream  function 
relaxation  factor 


for  the  case  of  pure  natural  convection 

quantities  at  {  —  A£ 

for  a  flat  plate 

conditions  at  the  wall 

condition  at  the  free  stream 


a 


■'M 
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INTRODUCTION 


Heat  transfer  by  natural  convection  along  a  vertical  cylinder  has  been  analyzed  rather 
extensively.  Elenbaas  (1948)  used  Iangmuir's  stagnant  film  model  to  evaluate  the  heat  transfer 
coefficients  for  vertical  cylinders  with  uniform  wall  temperature  (UWT).  Sparrow  and  Gregg 
(1956)  reworked  this  same  problem  using  a  power  series  expansion  method  for  Prandtl  numbers, 
Pr,  of  0.72  and  1  for  curvature  parameter,  AN  =  2-£-Gr;l/4,  up  to  about  1.5.  Kuiken  (1968)  also 

*o 

used  a  power  series  solution  and  investigated  a  power  law  variation  in  the  wall  temperature  and 
the  case  of  uniform  surface  heat  flux  (UHF)  for  0.7  <;  Pr  ^  100.  Ihc  power  series  solution  was 
also  applied  by  Fujii  and  Uehara  (1970)  to  cylinders  with  variable  wall  temperatures  for  AN  up 
to  2.73  and  for  0.72  ^  Pr<  100.  To  reduce  truncation  errors  of  the  power  series  and  the  local 
similarity  solutions,  Minkowycz  and  Sparrow  (1974)  employed  the  local  nonsimilarity  method 
and  obtained  results  for  Pr=  0.733  and  0  <  AN  <  10. 

Recently,  Lee  et  al.  (1986a,  1987)  analyzed  natural  convection  along  slender  vertical 

cylinders  under  UWT  and  UHF  conditions  as  part  of  their  mixed  convection  study  by 

employing  a  weighted  finite-difference  method  designed  to  overcome  the  numerical  difficulties 

associated  with  increasing  surface  curvature.  They  presented  results  for  0.1  <  Pr<  100  that 

cover  An  up  to  50  for  the  UWT  case  and  A*N  =  2-£-Gt*~1/s  up  to  50  for  the  UHF  case.  More 

*0 

recently,  Lee  et  al.  (1988)  analyzed  natural  convection  along  a  vertical  cylinder  for  the  case  of 
power  law  variation  of  wall  temperature,  Tw(x)  =  T^  +  ax",  and  found  that  the  results  of  Lee  et 
al.  (1986a)  were  not  accurate  for  AN  >  10  because  of  the  improper  step  sizes  used  in  the  radial 
direction.  They  presented  results  for  0. 1  <  Pr  <  100,  covering  AN  up  to  50  and  n  =  -0.4,  -0.2,  0, 
0.2,  and  0.5.  Lee  et  al.  (1987)  attached  a  qualifying  statement  to  their  study  of  the  UHF  case 
as  to  the  accuracy  of  their  results  for  curvature  parameters  A*N  >  10  and  attributed  this  also  to 
the  use  of  improper  step  sizes  in  rj.  It  appears  that  no  experimental  results  for  natural 
convection  along  vertical  cylinders  under  nonuniform  heating  have  been  reported  in  the 


literature. 


To  date  no  analytical  study  seems  to  have  been  presented  for  the  case  of  variable  surface 
heat  flux  along  a  vertical  slender  cylinder  and  this  has  prompted  the  present  study. 
Furthermore,  since  the  results  of  Lee  et  al.  (1987)  are  not  correct  for  large  values  of  A*N,  they 
need  to  be  reassessed. 


ANALYSIS 

Consider  a  semi-infinite,  vertical  cylinder  with  radius  r„  that  is  aligned  in  a  quiescent 
ambient  fluid  at  temperature  TM.  The  axial  coordinate  x  is  measured  upward  for  qw  >  0  and 
downward  for  q„  <  0,  while  the  radial  coordinate  r  is  measured  from  the  axis  of  the  cylinder. 
The  surface  of  the  cylinder  is  subjected  to  an  arbitrary  heat  flux  qw(x),  and  the  gravitational 
acceleration  g  is  acting  downward  The  fluid  properties  are  assumed  to  be  constant  except  for 
variations  in  density  which  induce  the  buoyancy  force.  By  employing  the  laminar  boundary 
layer  assumptions  and  making  use  of  the  Boussinesq  approximation  the  governing  conservation 
equations  for  the  problem  under  study  can  be  written  as: 

Continuity. 

~  (m)  +  z. — (rv)  =  0  ( I ) 


Momentum: 


+  S0(T-Too) 


Fnergy: 


dx 


+  v 


dT 

dr 


(2) 


(3) 


In  these  equations  u  and  v  are  the  velocity  components  in  the  x  and  r  directions 


respectively;  T  is  the  fluid  temperature;  and  v,  f}  and  a  are,  respectively,  the  kinematic  viscosity, 
the  volumetric  coefficient  of  thermal  expansion,  and  the  thermal  diflusivity  of  the  fluid. 


The  boundary  conditions  are 


u  =  v  =  d, 

OT 

-  qwM 

k  -r-'o 

(4a) 

u  — ♦  o ,  r  -►  i 

OC  as  r  -  «> 

(4ft) 

u  =  0,  1  =  1'^ 

at  x  =  0,  r  >  r0 

(4c) 

In  writing  the  last  boundary  condition  (4c)  it  is  assumed  that  the  flow  and  thermal  boundary 
layer  thicknesses  are  zero  at  the  leading  edge  of  the  surface. 


To  proceed  with  the  analysis,  the  conservation  equations,  along  with  the  boundary 
conditions,  are  transformed  into  a  dimensionless  form  by  introducing  the  following 
dimensionless  variables: 


1  -^(ttr//5),/5.  X  =  -f2L(C5rx*/5)- 

-r„x  *o 


f(T  r,)  =  *(xfr)/[5vr0(Gr//5),/5],  0(X.  rj)  = 


(  I  _  T  )(Grx*/5) 


qw(*)  */k 


where  Gr„*  =  ,g/?qw  (x)x4/AvJ  is  the  modified  Grashof  number,  >/  is  the  pseudo-simi!  irity  variable, 
f(/l,  rj)  is  the  reduced  stream  function,  0(X,  rj)  is  the  dimensionless  temperature,  X  is  the  curvature 
parameter,  and  i/<(x,  r)  is  the  stream  function  that  satisfies  the  continuity  equation,  with 
u  =  (<  i/>/dr)/r  and  v  =  —  (rty/dx)/r. 


The  transfonnation  yields 

(I  f  >,Z)f'"  4  Z.f'  4  (y  4  4)fP'  -  (2y  +  3)P2  +  0  -  A(y  -  1)(V'-^-  -  P— ) 

\  <7/(  OA  f 


;i  v  ,,X)U"  4  XO'  4  I’r(y  +  4|W'  -  Pr(4y  4  l)ffl  l»r/l(y  -  1  )(<)'  '1f  -  P^-\ 

V  d/  dyi  / 


f(Z,0)  -  P(T  ())-(),  ff'(i,0)  =  -  1 

P(T  tx>)  —  o,  d(  z ,  uo)  —  o 
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q*  dx 


and  the  primes  denote  partial  differentiation  with  respect  to  >7. 

The  system  of  equations  (7)-(9)  along  with  (10)  represents  the  general  form  of  the 
transformed  boundary  layer  equations  under  arbitrary  surface  heat  flux,  qw(x).  For  the  case  of 
power  law  variation  qw(x)  =  ax",  one  has  from  equation  (10)  that 


For  long  slender  cylinders,  the  curvature  parameter  X  =  2(x/r0)(Gr//5)  1,5  can  be  large.  lo 
lower  the  maximum  value  of  calculations  in  the  x  or  X  coordinate,  one  introduces  a  new  £(x) 


.■ariablc  defined  by 


£  =  [(x/r0)Grx*~'/5],/2 


such  that  X  —  C  with  C  =  2  •  5I/5. 


Substituting  equations  (11)  and  (12)  into  equations  (7)-(9)  results  in: 


Momentum: 


(1  +  a „)P"  +  a ,f'  +  a2ff”  +  a3P2  +  a 40  =  a5^f'-^-  -  f-f-) 


Fnergy: 


(I  +  a,^)£7"  +  a  ,0'  +  Fr  a2f<?'  +  Pr  a  6P0  =  Fr  ^(o~  - 


Boundary  Conditions  : 


f(£,  0)  =  f'(s,  0)  =  0,  0'(£,<» 

r({,oo)  =  0(^,c»)  =  () 


where 


■M 
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13 


Nux  =  —  =  -=r—^ - - 

x  l  T  —  T  i 
*  1  w  1  oo  ^ 


which  can  be  expressed  in  the  form 


Nux(Gr//5)-I/s  =  -  *'({, 0) 


The  average  Nusselt  number  is  defined  by  NuL  =  hl./k  where  the  average  heat  transfer 


coefficient,  h,  is  obtained  from 


This  leads  to 


-±r 

iJn 


NuL(GrL*/5r1/5  =  - 0)d{ 

i  ~  n  jn 


where  {L  =  {  at  x  =  L,  M  =  (7  +  3n)/(l  -  n),  and  N  =  (8  +  2n)/(n  -  1) 


Next,  from  the  definition  of  local  wall  shear  stress  r,  =  j/(du/dr)r,  ,0  one  obtains 


rw=5^GrxV5)3/5r'(i,0) 

x 


'i  nc  axial  velocity  distribution  can  be  written  as 


=  5(Gr//5)2/5f  «,n) 


and  the  temperature  profile  is  given  by  equation  (17). 


METHOD  OF  SOLUTION 


Equations  ( 1 8)-(20)  constitute  a  system  of  nonlinear  partial  differential  equations  in  the 
(i,rj)  coordinates  with  parameters  Pr  and  n.  The  method  described  in  Lee  et  al.  (1986b)  is 


■'V-V 

v 

•  ,W>V 


■m**. 
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employed  to  solve  this  system  of  equations.  The  first  step  is  to  convert  the  terms  involving  d/dt, 


in  the  following  manner. 


<311  P  /ii  tj  v  /  <311  ^  P  /it  it  \ 

nr  -  nr ("“,y nr 


where  H  is  a  function  of  (<J,  rj),  H„  =  H0  +  (q/p)A£(dH/d£)0,  and  the  subscript  "o  denotes 
quantities  at  £  -  A£.  The  values  for  p  and  q  arc  p  =  1  and  q  Oat  {  =  0  and  ^  =  A^,  and  p  =  2 
and  q  =  1  thereafter  for  £  >  2A£. 

The  f  and  4>  equations  arc  then  linearized  by  using  the  simple  technique  for  the  product  of 
two  arbitrarv  functions  F  and  G  ,  defined  as 


FG  =  FG  +  GF  -  FG 


where  F  and  G  indicate  values  to  be  guessed.  The  quasi-lineariz.cd  function  FG  expressed  by 

equation  (29)  will  provide  quadratic  and  monotone  convergence.  By  employing  equation  (28) 

for  the  terms  and  <£'(,<;,  0)  •—(  »,  j  and  applying  (29)  to  the  terms  ff", 

°Z  o£\<t>  (£,  0)  / 

P2,  f<£',  and  f  <£,  one  obtains 

A0P"  4-  (3|  +  A[)P'  +  A2P  +  A3f  +  A4<£  =  As  (30) 


n0<fi"  +  (a,  +  B  ,)<*>'  +  B2<£  +  BjP  +  n4f  =  Bs 


where 


A0=l+aIr/l  A,  =  a*f  +  d*f0  A2  =  2b  *f'-d*P0 
A3  =  a*P',  A4  =  a7  ,  A5  =  a*fP'  +  b*f '2 

B0  =  Aq  ,  B,  =  Pr  A,  ,  B2  =  Pr  e*f'.  B3  =  Pr  (c**  -  d*tQ) 
B4  =  PraV',  Bs  =  Pr(a*f<5>'  +  e*f'<J>) 


a*=a2-d*,  b*  =  a3+d*,  d*  =  a5  p/A<j;  ,  c*  =  a6  +  a8d* 
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and  the  other  variables  are  as  previously  defined. 

The  system  of  equations  has  now  been  reduced  to  a  system  of  quasi-linear  ordinary 
differential  equations.  By  defining  g  =  f  and  applying  the  weighting  factors  (see  Lee  et  al., 
1986b),  one  can  write  equations  (30)  and  (31)  in  the  finite-difference  form  as: 


f,  -  f(-_,  -  (A,,/2)(g,.  +  g(_,)  =  0 


(36) 


Equations  (36 )-(40)  constitute  a  system  of  algebraic  equations  that  can  be  written  in  the 
matrix  form  as 

[A][X]  =  [B]  (41) 

where  [A]  is  a  band  matrix  of  order  3n  and  bandwidth  seven.  The  array  [X]  which  contains 
the  solution  in  the  form  (f,,  g,,  <$>u  f2,  gj,  <t>2,  ...,  fn,  g„,  4>n)T  is  a  column  matrix  of  order  3n. 
The  matrix  [B]  is  a  column  matrix  of  order  3n  which  contains  the  right  hand  sides  of  equations 
(36)-(38)  and  (40).  The  matrix  [A]  is  approximately  diagonally  dominant  and  equation  (41)  can 
be  solved  by  the  Gaussian  elimination  technique  with  high  accuracy. 

To  obtain  values  for  f"  and  <f>',  the  results  for  f  and  <f>  along  with  the  boundary 
conditions 

n(,0)  =  -a1r(U)-a4,  <*>"(£,  0)  =  -  a,*'(£,0)t  P"(6  <»)-*"«.  oo)  -  0  (42) 

are  used  in  a  cubic  spline  interpolation  routine  (see,  for  example,  Burden  and  Faires  (1985)).  As 
the  solution  converges,  the  boundary  conditions  become  exact  and  the  values  for  f '  and  <f>'  can 
be  obtained  with  high  accuracy. 

Although  equation  (29)  should  provide  quadratic  monotone  convergence,  difficulties  arose 
as  the  surface  curvature  was  increased.  This  is  due  to  the  transformation  to  the  <f>  variable  with 
the  resulting  term  in  equation  (35).  Here  the  guess  of  <£'(£,  0)  is  needed  and  for  high  values  of 
surface  curvature  simply  updating  <f>'(£,  0)  as  the  next  guess  may  not  be  sufficient  to  obtain 
convergence  of  solutions.  To  improve  convergence,  attempts  were  made  to  improve  the  initial 
guesses  of  f,  f ,  f",  <f>,  and  <f>'  at  each  £  value. 

The  present  method  of  solution  employs  a  quasi-linearization  of  the  original  nonlinear 
system  of  equations  and  requires  initial  guesses  for  f,  f,  f",  <f>,  and  <p'.  The  flat  plate  solution 
(i.c.,  £  =  0)  for  the  uniform  surface  heat  flux  case  (UHF)  was  obtained  for  I’r  =  0.7  and  these 
results  were  used  as  the  initial  guesses  for  all  other  combinations  of  Pr  and  n  at  £  =  0.  For 
£  =  A£,  initial  guesses  were  taken  to  be  the  results  at  £  =  0.  At  £  =  2A£,  a  linear  extrapolation 
of  the  results  at  £  =  0  and  £  =  A£  was  used.  For  £  i>  3A£  a  three  point  inverse  polynomial 
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extrapolation  was  used  (see,  for  example,  Traub  (1964)).  This  approach  was  found  to  improve 
convergence  of  solutions  up  to  40%  faster  than  simply  using  the  previous  node's  values  as  the 
first  guess. 

To  improve  convergence  of  solutions,  it  was  found  that  the  value  needed  to  be  taken 
larger  as  the  value  of  the  surface  curvature  parameter  £  was  increased.  The  value  of  ^was 
initially  set  to  10  for  Pr  >  0.7  and  15  for  Pr  =  0.1.  It  was  increased  by  5  when  the  value  of 
P  or  <f>,  at  r\  =  O^Sf/^was  greater  than  0.001.  Since  P(£,  oo)  =  <f>(£,  oo)  =  0,  the  P  and  <j>  values 
greater  than  0.001  at  r\  =  0. 98>j00 indicate  that  the  corresponding  boundary  layer  thickness  extends 
past  To  make  this  adjustment,  all  values  of  f  and  <j>  and  their  derivatives  were  set  to  zero  for 
Tl>rl<x,  except  the  values  for  f  were  taken  to  be  the  value  at  the  old  extended  out  to  the  new 
This  adjustment  introduced  less  than  1%  difference  in  the  values  of  the  local  Nusselt 
number  when  compared  to  a  constant  run  at  rj^  =  30.  This  error  was  due  primarily  to 
initially  being  less  than  30  rather  than  to  the  adjustment  itself. 

The  solution  method  is  an  iterative  scheme  and  a  solution  was  considered  to  be 
convergent  when  the  calculated  values  for  f,  P,  and  <j>  differed  from  the  last  guess  of  the 
respective  values  by  less  than  10  4  at  all  nodes  (i.e.,  at  all  r\  values  for  a  given  £).  When  these 
criteria  failed,  new  guesses  for  f,  P,  and  <f>  were  found  using  a  weighted  average  of  the  last  guess 
and  the  resulting  calculation.  That  is 


^new  —  «f+(l  ^jfold’  ^  new  —  wP  +  (l  ^jf  old,  ^new  —  w  $  T  ( 1  “O'f’old  (43) 


in  which  to  is  a  relaxation  factor.  Generally  w  =  1  resulted  in  quick  convergence.  However,  if 
r\  «,was  too  small  or  for  some  small  values  of  Pr  and  n  it  was  sometimes  necessary  to  use 
under-relaxation  and  set  w  =  0.5  to  facilitate  convergence. 

It  was  found  that  as  £  was  increased,  errors  resulted  due  to  an  increase  in  the  boundary 
layer  thicknesses.  By  using  a  step  size  of  Ar/  =  0.01  errors  were  reduced  for  P'(£,  0)  and 
<f>'(£,  0)  at  high  values  of  the  curvature  parameter  £.  It  was  also  found  that  the  solution  was  not 
sensitive  to  the  step  size  in  £  and  A£  =  0. 1  was  used. 


V.jfc 


RESULTS  AND  DISCUSSION 


Typical  numerical  results  for  the  case  of  power  law  variation  of  surface  heat  flux  q„(x) 
were  obtained  for  values  of  the  exponent  n  of  -0.5,  -0.25,  0,  0.25,  and  0.5  for  Prandtl  numbers  of 
0.1,  0.7,  7,  and  100  over  the  range  of  0  ^  {  <  5  (i.e.,  0  <  A*N  <  50).  Representative  values  of  n 
fall  within  the  physical  limits  of  realism,  —  1  <  n  <  1,  for  power  law  variation  of  the  surface  heat 
flux  [see,  for  example,  Gebhart  (1971)]. 

lee  et  al.  (1987)  qualified  their  results  at  high  values  of  the  surface  curvature  as  a  result  of 
using  an  improper  step  size  Arj,  so  it  is  worth  noting  how  step  size  affects  the  local  Nusselt 
number.  For  the  UHF  case  with  A*N  =  50  and  Pr  =  0.7,  the  present  calculation  yields 
Nu,Gr*~I/J  =  16.6  using  Arj  =  0.036  and  =  1.24.  This  compares  to  NuIGr*~l/5  =  6.43  for 
At]  =  0.01  and  =  45.  This  trend  is  similar  to  what  Lee  et  al.(  1988)  found  for  the  UWT  case. 
It  was  found  that  decreasing  the  step  size  or  increasing  lowered  the  local  Nusselt  number. 
This  finding  explains  why  valid  results  were  obtained  at  small  values  of  the  surface  curvature 
where  the  range  of  r\  is  small  and  errors  appeared  for  large  values  of  curvature  where  the  range 
of  ij  has  increased.  When  that  happens  and  the  boundary  conditions  are  imposed  at  a  point 
that  is  actually  within  the  boundary  layer,  the  flow  and  thermal  profiles  are  lowered  thereby 
decreasing  {”(£,  0)  and  <£'(£,  0).  As  0)  becomes  more  negative  the  local  Nusselt  number 
increases. 

With  increasing  for  increasing  surface  curvature  it  becomes  necessary  to  determine  the 
degree  to  which  the  boundary  layer  assumptions  remain  valid.  For  boundary  layer  theory  to  be 
valid  (rt  -r„)/x  is  of  the  order  of  magnitude  of  the  term  (Gr//5)_,/5  where  (r4  —  r„)  =  S  is  the 
thickness  of  the  boundary  layer.  The  equations  for  2  and  *\  can  be  combined  to  give 

Y  =  -i<Gr//5)’/5  =  i-[(I  +  ^)'/2  _  ,j  (44) 

where  Y  is  a  measure  of  the  boundary  layer  thickness. 

Typically  i\t  is  defined  as  the  value  of »;  at  r  =  r4(  =  r0  +  5)  where  f  and  <f>  are  within  0. 1  % 
of  their  maximum  value  within  the  boundary  layer.  For  X  =  1  this  is  typically  about  tj(  =  20 
and  the  right-hand  side  of  equation  (44)  is  about  7.  If  [(r4  -  r0)/x](Gr//5),/s  remains  constant  at 
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7  and  2  increases  to  50,  rjt  will  increase  to  about  620.  This  shows  that  the  range  of  ^  can 
extend  very  large  as  the  surface  curvature  increases  and  the  boundary  layer  assumptions  will 
remain  valid.  It  also  verifies  the  analysis  presented  in  I^ee  et  al.  (1988)  and  Minkowyc/.  and 
Sparrow  (1974)  for  free  convection  under  UWT,  that  the  term  that  is  a  measure  of  the  boundary 
layer  thickness  may  actually  decrease  despite  an  increase  in  tj6.  This  would  happen  in  the 
example  if  tj,  is  less  than  620  but  greater  than  20. 

An  interesting  finding  is  that  a  decrease  in  the  step  si/e  Arj  (i.e.,  an  increase  in  the  number 
of  nodes)  has  a  large  influence  on  reducing  errors  in  the  local  Nusselt  number  at  large  values  of 
the  surface  curvature.  Ibis  is  due  to  the  effect  of  spreading  the  error  imposed  at  the  edge  of  the 
boundary  layer  among  a  greater  number  of  nodes,  thereby  lessening  its  effect  on  the  values  of  <f>' 
at  the  wall.  In  actuality  then,  the  error  of  I/;e  ct  al.  (1987)  arose  from  the  use  of  that  was 
not  sufficiently  large.  Since  a  decrease  in  the  step  si/e  has  an  effect  on  the  accuracy  of  the  local 
Nusselt  number  at  high  surface  curvatures,  a  small  Atj  was  used  in  conjunction  with  an  increase 
in  as  was  explained  in  the  method  of  solution. 

Numerical  results  from  the  present  analysis  will  now  be  presented.  Attention  is  first 
turned  to  the  local  Nusselt  number.  For  given  values  of  Pr  and  n,  as  the  surface  curvature  is 
increased  the  value  of  NuIGr*„  ,/5  is  found  to  increase.  Ibis  can  be  seen  in  Fig.  1  which  is  a  plot 
of  (Nux/Nu,J1)UHF  vs.  A*n  for  the  UIIF  case,  where  Nu^,  is  the  local  Nusselt  number  for  a 
vertical  flat  plate.  Ibis  plot  demonstrates  that  the  slope  of  such  a  normalized  plot  is  greater  for 
lower  Prandtl  numbers;  that  is,  the  relative  effect  of  surface  curvature  on  the  Nusselt  number  is 
stronger  for  fluids  with  a  smaller  Prandtl  number.  It  also  shows  a  nearly  linear  relationship 
which  is  useful  for  correlations  presented  later.  Table  1  presents  the  values  used  in  plotting  Fig. 
1  prior  to  normalization.  It  should  be  noted  from  the  table  that  as  the  curvature  increases  the 
effect  of  Prandtl  number  on  the  Nusselt  number  diminishes.  Ibis  can  be  explained  by  noting 
that  as  the  curvature  increases  the  first  two  terms  of  the  energy  equation  becoir  dominant  and 
the  problem  becomes  essentially  independent  of  Pr. 
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Figure  2  is  a  plot  of  (Nu,/Nuum  )p  vs.  n.  demonstrating  the  effeet  of  n  on  the  local  Nusselt 
number  for  the  vertical  flat  plate.  It  shows  that  a  lower  value  of  n  results  in  a  lower  value  of 
Nu,Gr*“1/s.  It  also  shows  that  for  an  increase  in  the  Prandtl  number  there  is  a  decrease  in  the 
deviation  of  the  Nusselt  number  from  the  L'HF  case. 

Table  2  shows  the  ratio  of  NuJN'u,  UHF  for  various  values  of  A*N.  It  can  be  seen  from  the 
table  that  as  the  curvature  increases  the  effect  of  n  on  Nu,  decreases  and  the  Nusselt  number 
ratio  approaches  1.0.  This  can  be  explained  again  by  the  dominance  of  the  first  two  terms  in 
the  governing  equations  for  large  values  of  the  curvature,  which  essentially  makes  the  system  of 
equations  independent  of  n.  Table  3  lists  the  average  Nusselt  number  results  NuLGr  V/5-  For 
given  values  of  Pr  and  n,  the  average  Nusselt  number  is  seen  to  increase  as  the  surface  curvature 
increases. 

The  P'((J,  0)  results  are  shown  in  Table  4.  It  can  be  seen  that  the  value  of  P'(£,0) 
decreases  as  the  surface  curvature  parameter  increases.  However,  a  decrease  in  P'(£,  0)  does  not 
necessarily  mean  that  the  local  wall  shear  stress  t„  decreases  as  x  increases  along  a  constant 
radius  cylinder.  Equation  [26]  can  be  rewritten  as 

rw  oc  0),  - 1  <  n  <  1  (45) 

For  n  2: 0,  the  power  to  which  {  is  raised  will  cause  tw  to  increase  despite  a  decrease  in 

P'({,  0).  As  n  decreases,  the  exponent  of  {  in  equation  (45)  decreases  until  at  n  =  -  2/3  it  will 
be  0  and  the  {  term  will  have  no  effect  on  t„  Since  P'(£,  0)  decreases  with  an  increase  in  the 
surface  curvature,  there  will  be  some  value  of  n  in  the  range  -  2/3  <  n  <  0  where  tw  will  decrease 
with  increasing  x  for  a  constant  radius  cylinder,  Therefore,  when  n  is  negative  it  is  possible  for 
rw  to  decrease  with  increasing  x.  This  is  particularly  true  for  low  Prandtl  numbers  for  which 
P'(£,  0)  decreases  at  a  greater  rate. 

The  observation  that  P'({,  0)  decreases  with  increasing  {  agrees  with  results  reported  by 

Kuiken  (1968).  It  is  interesting  to  note  that  this  trend  is  opposite  to  that  for  pure  forced 

convection  or  for  pure  natural  convection  under  the  power  law  variation  of  temperature.  To 
understand  why  this  occurs  physically,  one  should  consider  an  increase  in  the  surface  curvature 
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parameter  to  imply  a  decrease  in  the  radius  of  the  cylinder  at  a  given  x  location.  For  constant  x, 
Tw(x)  is  constant,  and  the  driving  force  for  pure  free  convection  under  the  variable  wall 
temperature  condition  is  constant.  For  constant  x  and  u„  the  driving  force  for  pure  forced 
convection  is  constant.  In  these  cases,  as  r0  decreases  the  velocity  near  the  wall  increases  and 
f ' '({,  0)  increases.  However  for  the  variable  heat  flux  case,  at  a  constant  x  location,  qw(x)  is 
constant.  When  qw  is  constant,  as  r0  decreases,  T„  -  TB  must  decrease  which  will  in  turn 
decrease  the  magnitude  of  the  velocity  profile  and  P'(£,  0).  This  explains  why  P'(iJ,  0)  increases 
for  the  cases  of  pure  forced  convection  and  for  natural  convection  under  the  power  law  variation 
of  surface  temperature,  but  decreases  for  natural  convection  under  the  power  law  variation  of 
surface  heat  flux. 

It  is  recognized  that  there  are  some  irregularities  in  the  values  of  P'(£,  0)  in  Table  4  that 
occur  for  higher  Prandtl  numbers.  This  is  due  to  a  failure  to  increase  >jtjo  appropriately  in  the 
numerical  computations.  Although  this  caused  little  effect  in  the  local  Nusselt  number  results 
due  to  the  small  step  size  used,  it  was  found  that  the  velocity  profile,  and  hence  the  value  of 
P'(f ,  0),  were  more  sensitive  to  the  choice  of  This  is  only  logical  since  f  starts  at  zero  at  the 
wall,  increases,  and  then  returns  to  zero  in  the  free  stream.  If  is  chosen  to  be  too  small,  the 
maximum  velocity  within  the  boundary  layer  will  be  lower  than  it  should  be.  Because  the 
momentum  boundary  layer  is  thicker  than  the  thermal  boundary  layer  for  high  Pr,  errors  in 
P'(£,  0)  will  result  if  is  large  enough  for  the  latter,  but  too  small  for  the  former.  This  does 
point  out  the  need  to  revise  the  criterion  for  increasing  ^  to  properly  account  for  the  boundary 
layer  thickness  in  numerical  calculations. 

For  practical  applications,  the  local  and  average  Nusselt  number  results  from  the  present 
calculations  in  the  ranges  of  0^A*N<50,  — 0.5  <  n  ^  0.5,  and  0.1  ^  Pr  <  100  can  be 
correlated  in  the  following  form 

NuxGr*~1/S  =  a(Pr)[A(A)  +  f,(Pr)A](l  +  VW)  (46) 


where 


a(Pr)  =  Pr2/5(4  +  9Pr1/2  +  10Pr)_,/5 


Xv> 


A(A)  =  1  4-  0.09A1/2 


f,(Pr)  =  (0.032  4-  0.176PrOJ84) 


V  =  {[0.328  4-  0.343exp(  -2.12Prl/5)]  -  0.195n}n 


W  =  exp[  -  (0.0265  4-  0.0907Pr-a444)A°'8] 


NuLGr*Ll/S  =  — «(Pr)[B(A)  4-  f2(Pr)A](l  +  V) 


where 


B(A)  =  1  +0.08A1/2 


f2(Pr)  =  (0.026  4-  0.14Pr-0'89) 


V  =  (4vw  —  nw)/(4  4-  nw) 


W  =  exp(  —  0.5A  ) 


For  notational  convenience,  A  in  the  above  equations  stands  for  A  =  A*N  =  2-p-Gr*, ,li.  In 


equations  (52)  -  (56)  A  is  calculated  at  x=  L. 


It  is  interesting  to  note  that  for  the  UHF  case  the  terms  V  and  V  in  equations  (46)  and 


(52)  become  zero  and  drop  out  of  the  equations.  It  is  also  interesting  to  note  that  for  the  flat 


plate  solution,  i.e.,  A  =  0,  the  terms  [A(A)  4-  f,(Pr)A]  and  [B(A)  4-  f3(Pr)A]  both  become  one. 


Therefore,  for  the  flat  plate  solution  under  UHF,  NuxGr'*xl/5  =  a(Pr),  where  a(Pr)  is  taken  from 


Fujii  and  Fujii  (1976).  The  maximum  error  in  the  correlations  for  the  local  and  average  Nusselt 


numbers  is  less  than  5%  for  the  UHF  case  and  less  than  8.3%  for  the  variable  heat  flux 


solution. 


s 


CONCLUSION 


In  this  paper,  the  problem  of  natural  convection  in  laminar  boundary  layer  flow  along 
slender  vertical  cylinders  has  been  investigated.  A  finite  difference  method  in  conjunction  with  a 
cubic  spline  interpolation  scheme  is  used  in  the  numerical  calculations  to  remove  the  difficulties 
associated  with  the  "stiffness"  of  the  governing  equations  for  high  values  of  the  curvature 
parameter.  Additionally,  a  smaller  radial  step  size  and  a  larger  value  of  were  used  to  correct 
reported  errors  in  the  Nusselt  number  at  high  values  of  the  curvature  parameter.  'Hie  major 
findings  from  these  results  can  be  summarized  as  follows: 

(1)  Errors  in  the  thermal  and  velocity  profiles  result  if  the  magnitude  of  is  insufficient. 
The  calculated  values  of  the  local  Nusselt  number  will  be  too  large  and  those  of  f”(£,  0)  will  be 
too  small.  These  errors  can  be  reduced  by  decreasing  the  radial  step  size  A rj  or  by  increasing 
»/„.  A  decrease  in  the  radial  step  size  is  more  effective  at  reducing  errors  for  the  local  Nusselt 
number  but  it  does  not  correct  the  entire  velocity  and  thermal  profiles. 

(2)  For  the  power  law  variation  in  surface  heat  flux,  the  local  surface  heat  transfer  rate 
increases  with  increasing  value  of  the  exponent  n.  It  also  increases  with  increasing  curvature  and 
increasing  value  of  Pr.  However,  as  the  curvature  increases  the  effects  of  n  and  Prandtl  number 
on  the  Nusselt  number  diminish. 

(3)  The  behavior  of  the  average  Nusselt  number  is  similar  to  that  for  the  local  Nusselt 
number  for  UHF.  TTxe  effect  of  n  and  Pr  on  the  average  Nusselt  number  decreases  as  the 
curvature  increases. 

(4)  The  value  of  f”(£,  0)  decreases  with  increr..  ng  curvature,  decreases  with  increasing  n, 
and  decreases  with  increasing  Prandtl  number.  The  local  wall  shear  stress  rw,  however,  increases 
as  x  increases  for  a  constant  radius  cylinder  with  n  >  0.  It  may  decrease  with  increasing 
curvature  for  n  <  0. 
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MIXED  CONVECTION  ALONG  A  VERTICAL  CYLINDER 
WITH  VARIABLE  SURFACE  TEMPERATURE 


J.  J.  Ilcckel,  T.  S.  ('hen,  and  B.  I:.  Armaly 
Department  of  Mechanical  and  Aerospace  Engineering 
University  of  Missouri -Rolla,  Rolla,  MO  65401 


ABSTRACT 


Mixed  convection  in  laminar  boundary  layer  flow  along  slender  vertical  cylinders  is 
analyzed  for  the  situation  in  which  the  surface  temperature  Tw(x)  varies  arbitrarily  with  the  axial 
coordinate  x.  It  covers  the  entire  mixed  convection  regime  from  pure  free  convection  (y  =  0)  to 
pure  forced  convection  (x  =  1),  where  y  =  [  1  +  (Grx/RejU,/4]  1  is  the  mixed  convection 
parameter.  The  governing  boundary  layer  equations  along  with  the  boundary  conditions  are 
first  cast  into  a  dimensionless  form  by  a  nonsimilar  transformation  and  the  resulting  system  of 
equations  is  then  solved  by  a  weighted  finite-difference  method  of  solution  in  conjunction  with 
cubic  spline  interpolation.  Sample  calculations  were  performed  for  the  case  of  power  law 
variation  in  surface  temperature,  Tw(x)  —  T„  =  ax",  for  fluids  with  Prandtl  numbers  of  0.1,  0.7, 
7,  and  100  over  a  wide  range  of  values  for  the  surface  curvature  parameter  0<;  A<  50  (or 
0  <  £  <,  5).  Local  and  average  Nusselt  numbers  as  well  as  the  local  wall  shear  stress  are 
presented.  It  is  found  that  the  local  Nusselt  number  in  the  form  Nu^ReJ/2  +  GrJ/4)  increases 
with  increasing  surface  curvature,  Prandtl  number,  and  the  exponent  n.  For  low  values  of 
constant  curvature  A,  the  local  Nusselt  number  initially  decreases  and  then  increases  as  y  goes 
from  0  to  1.  As  curvature  increases  a  linear  relationship  is  found  to  exist  between  the  Nusselt 
number  and  the  mixed  convection  parameter,  particularly  for  low  Prandtl  numbers.  The 
velocity  gradient  at  the  wall  is  found  to  increase  with  increasing  curvature,  decreasing  Prandtl 
number,  and  decreasing  n.  Correlation  equations  for  the  local  and  average  Nusselt  numbers  are 
also  presented. 
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NOMHNCLATURK 

local  skin  friction  coefficient,  2rjpu2x 

average  skin  friction  coefficient,  (2/pu^L) J^rwdx 

reduced  stream  function,  (i/</vr0)(Rei/2  +  GrJ'4)-1 

df/dt)  or  gravitational  acceleration 

local  Grashof  number,  g/?[Tw(x)  -  T^x3/*2 

modified  local  Grashof  number,  g/?qw(x)x4/A:v2 

Grashof  numbers  defined,  respectively,  as 
g/?[Tw(L)  -  TJL3/v!  and  g/?[Tw(r0)  -  TJiJ/v’ 

local  heat  transfer  coefficient  ,  qw/(Tw  -  T„) 

average  heat  transfer  coefficient 

thermal  conductivity  of  the  fluid 

an  arbitrary  length  of  cylinder 

exponent  used  in  the  mixed  convection  correlation 

exponent  in  the  power  law  variation  of  the  surface  temperature 
or  surface  heat  flux 

local  Nusselt  number,  Ax/A: 

average  Nusselt  number,  A  L/A 

Prandtl  number,  v/a 

radial  coordinate 

radius  of  the  cylinder 

Reynolds  number  based  on  x,  uMx/v 

Reynolds  number  based  on  L  and  r0  ,  respectively,  as  u.J./v  and  u,jjv 

fluid  temperature 

axial  velocity  component 

radial  velocity  component 

axial  coordinate 

transverse  coordinate 

dimensionless  axial  coordinate,  x/r0 


thermal  diffusivity 


m  m 


I 


/? 

y 

9 

An 

A*n 

Af 

A 


t 

P 

T- 

X 

X * 

n* 

n% 

a 


Subscripts 

F 

N 

o 


volumetric  coefficient  of  thermal  expansion 
variable  temperature  parameter,  [x/(Tw  -  T„)][d(Tw  -  T^J/dx] 
pseudo-similarity  variable,  [(r2  -  r2)/2r0x](Rc;/2  +  GrJ/4) 
dimensionless  temperature,  (T  -  T00)/[Tw(x)  -  T^] 

curvature  parameter  for  natural  convection  with  variable  surface  temperature, 
(2x/r0)Grx  1/4 

modified  curvature  parameter  for  natural  convection  with  variable  surface 
heat  flux,  (2x/r0)Gr%*£ 


curvature  parameter  for  pure  forced  convection,  (2x/r0)Re„ 


1/2 


curvature  parameter  for  mixed  convection  with  variable  surface  temperature, 
(2x/r0)(ReJ/2  +  Gri/V 

curvature  parameter  for  mixed  convection  with  variable  surface  heat  flux, 
(2x/r0)(Rei/2  +  Gr4^)' 

dynamic  viscosity 

kinematic  viscosity 

dimensionless  axial  coordinate,  [A/(2y)]1/2  or  [x/(r0Re0)]l/4 

fluid  density 

local  wall  shear  stress 

stream  function 

mixed  convection  parameter  for  variable  surface  temperature,  ( 1  +  i2J/4)  1 
mixed  convection  parameter  for  variable  surface  heat  flux,  (1  +  1 

buoyancy  parameter  for  variable  surface  temperature,  Gr,/ReJ 
buoyancy  parameter  for  variable  surface  heat  flux,  Gr4,/Re$/2 
modified  buoyancy  parameter,  [GrJReJ,1  n)],/4 
relaxation  factor 


for  the  case  of  pure  forced  convection 
for  the  case  of  pure  natural  convection 
quantities  at  £  — 


INTRODUCTION 


There  have  been  relatively  few  studies  on  mixed  convection  along  a  vertical  cylinder. 

Chen  and  Mucoglu  [1,2]  were  the  first  to  analyze  the  effects  of  bouyancy  forces  on  forced 

convection  along  vertical  cylinders  for  the  case  of  uniform  wall  temperature  (UWT)  and  uniform 

surface  heat  flux  (UHF).  They  used  the  local  nonsimilarity  method  of  solution  to  obtain  heat 

transfer  results  that  cover  the  surface  curvature  parameter,  AF  =  2— Re„1/2,  from  0  to  8  for 

1*0 

Prandtl  numbers  of  0.7  and  7.  Bui  and  Cebeci  [3]  analyzed  the  UWT  case  for  Pr  =  0.1,  1.0, 

and  10  for  curvatures  up  to  AF  =  10  using  a  central  difference  finite-difference  method  of 

solution.  The  studies  by  Chen  and  Mucoglu  [1]  and  by  Bui  and  Cebeci  for  the  UWT  case  used 

a  mixed  convection  parameter  fl*  =  Grx/Re2  which  varies  from  zero  for  pure  forced  convection 

to  infinity  for  pure  free  convection.  Mucoglu  and  Chen  [2]  presented  the  UHF  case  with  a 

mixed  convection  parameter  =  Gr*JResJ2.  Lee  et  al.  [4]  treated  the  UWT  case  for  slender 

cylinders  using  a  new  mixed  convection  parameter  y  =  ( 1  +  f2]'4)  1  which  varies  from  zero  for 

pure  free  convection  to  one  for  pure  forced  convection.  In  Lee  et  al.  [5]  the  UHF  case  is 

presented  using  x*  —  (I  +  D*i/5)'‘  which  also  varies  from  zero  for  pure  free  convection  to  one 

for  pure  forced  convection.  The  main  difficulty  encountered  in  these  studies  has  been  the 

numerical  solutions  for  large  surface  curvatures.  The  results  of  Chen  and  Mucoglu  [1,2]  agree 

well  with  those  of  Lee  et  al.  [4,  5]  up  to  surface  curvatures  A  =  2— (Re]/2  +  GrJ/4)  1  for  UWT  or 

A*  =  2^-(Re[/2  +  Gr*’75)'1  for  UHF  of  8  for  Pr  =  0.7  and  7.0.  The  results  of  Bui  and  Cebeci 
*0 

are  about  36  percent  higher  for  large  values  of  the  surface  curvature  and  Lee  et  al.  [4]  attributed 
this  to  the  latter's  use  of  a  central  difference  finite-difference  method  of  solution  which  has 
difficulty  handling  the  boundary  layer  equations  as  they  become  "stiff"  with  increasing  curvature. 
Lee  et  al.  [4,  5]  advocates  the  use  of  a  weighted  difference  finite-difference  method  [6]  along 
with  a  cubic  spline  interpolation  method  [7]  to  overcome  these  difficulties  associated  with  high 
values  of  curvature  or  Prandtl  number. 

Another  advantage  of  the  method  used  by  I  ee  et  al.  [4,  5]  was  that  it  was  applied  to  the 
full  range  of  mixed  convection  to  include  pure  free  convection.  The  case  of  natural  convection 
along  a  vertical  cylinder  has  been  extensively  studied  [8-  13],  Llenbaas  [8]  used  l^uigmuir's 
stagnant  film  model  to  evaluate  the  heat  transfer  coefficients  for  vertical  cylinders  under  UWT 


conditons.  Sparrow  and  Gregg  [9]  reworked  this  same  problem  using  a  power  series  solution 
for  Pr  =  0.72  and  1  covering  surface  curvatures  (AN  =  2-£-Gt~1/<)  of  about  1.5.  Kuiken  [10] 
also  applied  the  power  series  expansion  method  to  investigate  some  nonuniform  wall 
temperature  conditions  and  obtained  results  for  0.7  <  Pr  <  10.  Similarly,  Fujii  and  Uehara  [11] 
treated  the  case  of  variable  wall  temperature  and  obtained  results,  again  by  the  power  series 
solution,  for  0.72  <  Pr  <,  100  covering  AN  from  0  to  2.73.  To  overcome  truncation  errors  of  the 
power  series  and  the  local  similarity  solution,  Minkowycz  and  Sparrow  [  1 2]  employed  the  local 
nonsimilarity  method  of  solution  and  obtained  results  for  Pr  =  0.733  covering  0  <  AN  <  10. 
Lee  et  al.  [4,  5]  examined  both  the  UWT  and  IJIIF  cases  of  natural  convection  along  a  vertical 
cylinder  as  part  of  their  mixed  convection  study  and  employed  a  modified  finite-difference 
method  designed  to  overcome  the  difficulties  associated  with  increasing  curvature.  As  a  result, 
they  were  able  to  obtain  results  up  to  AN=  50  for  0.1  <  Pr<  100.  Very  recently,  Lee  et  al. 
[13]  treated  the  variable  wall  temperature  problem  for  natural  convection  along  a  vertical 
cylinder.  They  found  that  the  results  of  [4]  were  inaccurate  for  AN  >  10  because  an  improper 
step  size  was  used  in  the  radial  direction  . 

The  significant  finding  for  pure  forced  or  pure  free  convection  under  a  variation  in  surface 
temperature  along  a  vertical  cylinder  is  that  as  the  surface  curvature  increases  at  a  given  x 
location  both  the  heat  transfer  coefficient  and  the  skin  friction  increase.  For  the  UHF  case  of 
natural  convection,  as  the  surface  curvature  at  a  given  x  location  increases,  the  heat  transfer 
coefficient  also  increases,  but  the  local  wall  shear  stress  decreases.  Minkowycz  and  Sparrow 
[12]  showed  that  the  error  caused  by  the  power  series  solution  increases  with  increasing 
curvature.  Results  by  lee  et  al.  [13]  showed  that  the  results  of  Fujii  and  Uehara  [11]  are  valid 
only  for  AN  up  to  about  1.5.  Lee  et  al.  [4]  demonstrated  that  the  central  difference  scheme  used 
by  Bui  and  Cebeci  [3]  produces  large  errors  as  the  curvature  increases,  but  the  results  of  [4,  5] 
for  A  >  10  are  not  accurate  because  they  used  too  large  a  step  size  in  the  radial  direction. 

For  comparison  purposes,  results  for  mixed  convection  along  a  vertical  flat  plate  are 

needed.  There  are  numerous  studies  on  mixed  convection  along  a  vertical  flat  plate  (see,  for 

example,  [14]  —  [19]  ).  The  Nusselt  number  results  for  this  flow  geometry  have  been  correlated 

m  m  m 

by  Churchill  [20]  using  the  simple  form  Nu  =  NuF  +  Nun  ,  with  m  =  3,  where  \’uF  is  the 


Nusselt  number  for  pure  forced  convection  and  Nun  is  the  Nusselt  number  for  pure  free 
convection.  Chen  et  al.  [14]  have  carried  out  extensive  correlations  for  mixed  convection  along 
vertical,  inclined,  and  horizontal  flat  plates  and  verified  them  with  experimental  data  for  air. 

From  the  past  studies,  it  is  clear  that  results  for  mixed  convection  under  UWT  conditions 
at  high  values  of  curvature  need  to  be  reassessed.  In  addition,  the  combined  effects  of  variable 
surface  temperature  and  curvature  on  the  flow  and  heat  transfer  for  mixed  convection  along 
slender  vertical  cylinders  have  not  been  studied  for  the  entire  mixed  convection  regime.  This  has 
motivated  the  present  study.  It  is  noted  that  free  convection  along  slender  vertical  cylinders  has 
been  analyzed  for  the  case  of  power  law  variation  in  the  wall  temperature.  Details  of  this 
analysis  can  be  found  in  [13]  and  for  completeness  its  highlights  are  also  given  in  Appendix  C. 
Some  results  for  this  problem  were  obtained  and  will  be  discussed  later.  An  analysis  of  the 
mixed  convection  problem  under  a  power  law  variation  in  the  wall  temperature  follows. 


ANALYSIS 


Consider  a  semi -infinite,  vertical  cylinder  with  radius  r0  that  is  aligned  parallel  to  a 
uniform,  laminar  free  stream  with  velocity  u^  and  temperature  T„.  The  axial  coordinate  x  is 
measured  in  the  direction  of  the  forced  flow  and  the  radial  coordinate  r  is  measured  from  the 
axis  of  the  cylinder.  The  surface  of  the  cylinder  is  subjected  to  an  arbitrary  variation  in 
temperature  Tw(x),  and  the  gravitational  acceleration  g  j  acting  downward.  Fluid  properties  are 
assumed  to  be  constant  except  for  variations  in  density  which  induce  the  buoyancy  force.  By 
employing  the  laminar  boundary  layer  assumptions  and  making  use  of  the  Boussinesq 
approximation,  the  governing  conservation  equations  can  be  written  as: 

Continuity: 


!_(ru)  +  (rv)  =  0 


(1) 


Momentum: 


Energy: 


(3) 


The  positive  sign  in  equation  (2)  applies  to  upward  forced  flow  and  the  negative  sign  to 
downward  forced  flow.  In  these  equations,  u  and  v  are  the  velocity  components  in  the  x  and  r 
directions,  respectively;  T  is  the  fluid  temperature;  and  v,  /?  and  oc  are,  respectively,  the  kinematic 
viscosity,  the  volumetric  coefficient  of  thermal  expansion,  and  the  thermal  diffusivity  of  the  fluid. 

The  boundary  conditions  are 

u  =  v  =  0  ,  T  =  T w(x)  at  r  =  r0  (4) 

u  -»  0  ,  T  -*  T^  As  r  ->  oo  (5) 

u  =  0  ,  T  =  T^  At  x  =  0,  r  >  r0  (6) 


In  writing  equation  (6)  it  is  assumed  that  the  flow  and  thermal  boundary  layer  thicknesses  are 
zero  at  the  leading  edge  of  the  cylinder  surface. 


The  conservation  equations  and  the  boundary  conditions,  are  then  transformed  into  a 
dimensionless  form  by  introducing  the  following  dimensionless  variables: 


n  = 


— - — (Re^2  +  Gr.1/4 


2r0x 


). 


(7) 


f(z,  r,)  =  ^{x,r)/[vr0(Rei/2  +  GrJ/4)],  0(z,  >,)  =  (T  -  T,J/[Tw(x)  -  TJ  (8) 


X  =  (1  +  «i/4)“ 


(9) 


where  Gr„  =  g^[Tw(x)  -  Tw]  x3/vJ  is  the  local  Giashof  number,  Rex  =  u^x/v  is  the  local 
Reynolds  number,  tj  is  the  pseudo-similarity  variable,  z  is  the  dimensionless  axial  coordinate, 
f(z,  »j)  is  the  reduced  stream  function,  0( z,  y)  is  the  dimensionless  temperature,  i^(x,  r)  is  the 


stream  function  that  satisfies  the  continuity  equation,  with  u  =  (<ty/dr)/r  and  v  =  -  (dipjd\)l r, 
is  the  buoyancy  parameter  which  varies  from  zero  for  pure  forced  convection  to  infinity  for  pure 
free  convection,  and  y  is  the  mixed  convection  parameter  which  varies  from  zero  for  pure  free 
convection  to  one  for  pure  forced  convection. 

The  transformation  yields: 

(l+,,A)P"  +  Af"  +  ^-(2  +  (l-y)(y4-  l))fT 

-  ~(1  -  x)(v  +  1)P2  ±  (I  -  xfo  -  -  /.(p'-f-  -  P-g-) 


(1  +  ,,A)0"  +  A  O'  4-  -^(2  +  (1  -  y)(y  +  l))fO' 


Pryf’0 


=  -  Pr 


f-r§) 


(II) 


f(z,  0)  =  f  (z,  0)  =  0,  0( z,  0)  =  1 

P(Z,  oo)  =  y2,  0( Z,  co)  =  0 


(12) 


where 


A  =  2-2L(R  c'J2  +  Gt'JY' 


(13) 


is  the  mixed  convection  curvature  parameter  and 


^-^-n'wW-Too] 


(14) 


The  primes  in  equations  ( 1 0)-(  1 2)  denote  partial  differentiation  with  respect  to  tj,  and  the  plus 
and  minus  signs  in  front  of  the  term  (1  —  y)40  in  equation  (10)  now  represent  buoyancy  assisting 
flow  and  buoyancy  opposing  flow,  respectively. 

The  system  of  equations  (10)-(12)  represents  the  genera]  form  of  the  transformed 


boundary  layer  equations  for  variable  wall  temperature  Tw(x)  along  vertical  cylinders  in  mixed 
convection.  For  the  case  of  power  law  wall  temperature  distribution,  Tw(x)  —  T„  4  a\",  one  has 
from  equation  (14)  that 


Equations  ( 10)-(  1 2)  contain  three  x-dcpcndent  parameters,  z(x),  A(x),  and  x( *)■  These 


parameters  can  be  related  to  a  single  x-dcpcndent  parameter  £(x)  defined  by 


(‘(it)'1*' 


A  =  2^z,  *  =  (1  +Q0£(l+n)rI 


where 


no  =  [Gr0/Re’"n]l/4 


with  Gr0  =  Gr,  and  Re0  =  Re,  for  x  =  r0  Also,  the  right-hand  sides  of  equations  (10)  and  (11) 


become,  respectively 


-I'H-rf)  M°t-e f) 


Rewriting  equations  ( 10)-(  1 2)  one  has 


Momentum: 


(1  4-  a^)f”  +  a,f'  +  a2lf'  +  a3P2  +  a 40  =  a5(p'-|-  -  P-|0 


Energy: 


(1  +  a,ij)0"  +  a,0'  +  Pr  a2f0'  +  Pr  a^PO  =  Pr  a5^'^-  -  r-fr) 


Boundary  Conditions  : 


f(£,  0)  =  P(£,  0)  =  0,  Otf,  0)  =  1 
P({,oo)=a7t  0(£,  oo)  =  0 


(22) 


a,=A,  a2  =  -i-(2  +  (l-*)(n  +  1))  ,  a3  =  ^{x  -  l)(n  +  1) 
a4  =  ±(l-z)4.  as  =  -{/4,  a6  =  -n,  a7  =  *2 


With  A  and  y  related  to  £,  the  functions  f  and  f?  in  equation  (20)-(22)  arc  functions  of 
({,»/)  and  depend  on  three  constant  parameters  n,  Pr,  and  £20.  These  equations  arc  now  in  the 
form  that  can  he  solved  by  the  method  proposed  by  Ice  et  al.  [f>]  See  Appendix  !■’  for  the 
solution  method  as  applied  to  equations  ( 20>-( 22).  l  or  the  purpose  of  comparisons,  the 
problem  of  mixed  convection  along  a  vertical  flat  plate  with  a  power  law  variation  in  the  wall 
temperature  is  also  solved.  The  governing  equations  and  boundary  conditions  have  exactly  the 
same  form  as  equations  (20)-(23)  with  the  exceptions  that  for  this  case 


>1  -  (|-)(Rei/2  +  Gri/4> 

i  =  Z,  A  =  0  ,  a5-(l  +n)(l— *)x/4 


It  should  be  noted  that  this  problem  must  be  solved  numerically  from  £  =  1  to  f  =  0. 

The  physical  quantities  of  interest  include  the  local  and  average  Nussclt  numbers,  the  local 
and  average  friction  factors,  the  axial  velocity  distribution,  and  the  temperature  profile.  The 
local  Nussclt  number  and  the  local  friction  factor  can  be  expressed  as 

Nux/(ReJ/2  +  Gri/4)  =  -  0'(£,  0)  (25) 


Cf^  =  — =  2z~3f'(4, 0)Rex 

PU  oo 


The  average  Nusselt  number  and  the  average  friction  factor  can  be  expressed  as 
NV(Re,'/2  +  Gr|/4)  =  -4*,  <^2  [ 

Ju  X 


CfLRei/2  =  8^2f5Lz“3P'(^,0)^ 

•/fl 


*1 


ww 


N  V 
a  *  »  *  * 

>>*Cv 
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where  fL  =  £  at  x  -  L,  and  xl  =  X  at  x=  L.  The  axial  velocity  distribution  can  be  written  as 


U  nt.n) 


and  the  temperature  profile  is  given  by  0(£,  rj)  =  (T  -  T00)/(TW  -  Tj. 


(29) 


METHOD  OF  SOLUTION 

Equations  (20)-(22)  constitute  a  system  of  coupled  nonlinear  partial  differential  equations 
in  the  (£,  rj)  coordinates  with  parameters  Pr,  n,  Q0.  The  solution  method  described  in  [6]  was 
employed  to  solve  this  system  of  equations. 

The  first  step  is  to  convert  the  terms  involving  dfdS  in  the  following  manner. 

f--^-(M-H„)-q(^-)0  =  ^-(ll-Ho)  HO) 

where  H  is  a  function  of  (£,  >j),  H0  =  H0  +  (q/p)A£(<5H/d£)0,  and  the  subscript  "o'  denotes 
quantities  at  £  —  A£.  The  values  for  p  and  q  are  p  =  1  and  q  =  0  at  {  =  0  and  £  =  A£,-  and  p  =  2 
and  q  =  1  thereafter  for  {  >  2A£.  The  f  and  0  equations  arc  then  linearized  by  using  the  simple 
technique  for  the  product  of  two  arbitrary  functions  F  and  G  ,  defined  as 

FG  =  FG  +  GF  -  FG  (31) 

where  F  and  G  indicate  the  values  to  be  guessed.  The  quasi-linearized  function  FG  expressed 
by  equation  (31)  will  provide  quadratic  and  monotone  convergence.  By  employing  equation 

af  art  an 

(30)  for  the  terms  — — ,  ——  and  ~  and  applying  equation  (31)  to  the  terms 

o S  os  di 

fP' ,  f  ■ 2 ,  f O'  and  f  0  one  obtains 

AqP*  +  (aj  +  Aj)f '  +  A2P  +  Ajf  +  A^P  =  A5  (22) 

Bfflff"  +  (al  +  B|)0'  +  »20  +  »3P  +  B4f  =  1*5  (33) 


where 
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Aq  =  1  +  a,»j ,  A,  =  a*f  +  d*f0>  A2  =  2b*f' -  d*P0 
A3  =  a*P',  A4  =  a4  ,  As  =  a*fP'  +  b*f'2 

B0  =  Ao  ,  B,  =  Pr  A,  ,  B2  =  Pr  e*f',  B3  =  Pr  (e*0  -  d*0o ) 
B4  =  Pr  a*0‘,  Bs  =  Pr(a*f0'  +  e*f‘ ’0) 


a*  =  a2  —  d*  ,  b*  =  a3  +  d*,  d*  =  a5p/A^,  e*  =  a6  +  d* 


The  coefficients  a,  through  a6  and  the  functions  f„,  f0,,  and  0o  are  as  previously  defined. 

The  f  and  0  equations  have  now  been  reduced  to  a  system  of  quasi-Iinear  ordinary 
differential  equations.  By  defining  g  =  P  and  applying  the  weighting  factors  (see  I,ee  et  al.  [6]), 
one  can  write  equations  (32)  and  (33)  in  the  finite  difference  form  as: 

0 -  ff_i  -  (A^/2)(g,  +  g/_,)  =  0  (37) 


Ao/(Al?)  (a_lfv_|  +  a0&  +  <*+}%/+))  +  A2 8/  +  Aj0  +  A4 Qi  -  A5 


Bo/(A>?)  +  PcPi  +  P+}&i+\)  +  +  B3g,  +  B4f,  -  B5  (39) 

where  A>j  is  the  step  size  in  the  tj  direction  and  the  subscript  "i"  refers  to  values  at  the  nodal 
point 

In  equations  (38)  and  (39)  the  weighting  factors  are  defmed  as: 

a_i  =  W((  —  z(_|/2) ,  “+i  =  WKzi+i/2)  >  ao  =  -“-i  -«+i 

/L,=W,(-z  V,/2),  /?+1=WKz  \+1/2),  Po^-P-i-P+i  (40) 

Wf(z)  =  z/[1  -exp(-z)],  z  =  Aijfa,  +  A^/Aq  ,  /.♦  =  A»j(a1  +  B,)/B0 

Equations  (37)-(40)  are  applied  to  the  interior  points.  On  the  boundaries,  equation  (22) 
applies  such  that 


“.  m •  *  "  *  J 


PM 

pm 
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f,=g!=0,  0,  =  1 


fn-fn-l  -— (gn  +  8n-l)  =  0n  =  O.  8n  =  a7 


where  the  subscript  n  is  the  number  of  nodes  in  the  rj  direction  and  a7  =  y2  for  mixed 


convection. 


Equations  (37)-(41)  constitute  a  system  of  algebraic  equations  that  can  be  written  in  the 


matrix  form 


[A][X]  =  [B] 


where  [A]  is  a  band  matrix  of  order  3n  and  bandwidth  seven.  The  array  [X]  which  contains 


the  solution  in  the  form  (f,,  g,,  0,.  f2,  g2,  02,  ...,  fn,  g„,  0„)T  is  a  column  matrix  of  order  3n. 


The  matrix  [B]  is  a  column  matrix  of  order  3n  which  contains  the  nght  hand  sides  of  equations 


(37)-(39)  and  (41).  The  matrix  [A]  is  approximately  diagonally  dominant  and  equation  (42)  can 


be  solved  by  the  Gaussian  elimination  technique  with  high  accuracy. 


To  obtain  values  for  f '  and  O',  the  results  for  P  and  0  along  with  the  boundary  conditions 


P"(£,  0)  =  -  a,P'(^,  0)  —  a4,  0'%  0)  =  -  a,0U  0),  P"(£,  «)  =  0"({, «)  =  0  (43) 


are  used  in  a  cubic  spline  interpolation  routine  (see,  for  example,  Burden  and  Faires  [7]  ).  As 


the  solution  converges,  the  boundary  conditions  become  exact  and  the  values  for  P'  and  O'  can 


be  obtained  with  high  accuracy. 


The  present  method  of  solution  employs  a  quasi-linearization  of  the  original  nonlinear 


system  of  equations  and  requires  initial  guesses  for  f,  P,  P',  0,  and  O'.  The  flat  plate  solution 


(i.e.,  £  =  0)  for  the  uniform  wall  temperature  case  (UWT)  for  Pr  =  0.7  was  obtained  and  these 


results  were  used  as  the  initial  guesses  for  all  other  combinations  of  Pr,  n,  and  Ll0  at  £  =  0.  For 


5  >  A£,  good  convergence  was  obtained  by  simply  letting  the  solution  at  the  previous  node  be 


the  guesses  for  the  next  node. 
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The  solution  method  is  an  iterative  scheme.  A  solution  was  considered  to  be  convergent 
when  the  calculated  values  for  f,  f ,  and  0  differed  from  the  last  guess  of  the  respective  values  by 
less  than  10  4  at  all  nodes.  When  these  criteria  failed  to  be  met,  new  guesses  for  f,  P,  and  0 
were  found  using  a  weighted  average  of  the  last  guess  and  the  resulting  calculation.  That  is 

fnew  =  w  f  +  ( 1  -  to)fold,  f 'new  =  w  P  +  ( 1  -  to)f  ‘ 'old,  0new  =  w  0  +  (1  -  a>)0old  (44) 


Generally  a>  =  1  resulted  in  quick  convergence.  However,  for  G0  =  2  it  was  often  necessary  to 
use  undcr-rclaxation  with  ui  =  0.5  to  facilitate  convergence. 


It  was  found  that  numerical  results  depended  upon  the  choice  of  jj^.  As  £  was  increased, 
it  was  necessary  to  increase  The  value  of  was  initially  set  to  1 5  and  was  increased  by  5 
when  the  value  of  |P  —  a7|  or  0  was  greater  than  0.001  at  y]  -  0.98^.  Since 
P(£,  oo)  —  a,  =  0(£,  oo)  =  0,  the  |  P  —  a7 1  and  0  values  greater  than  0.001  at  »j  =  0.98^  indicate 
that  the  corresponding  boundary  layer  thickness  extends  past  To  make  this  adjustment  for 
free  convection,  all  values  of  f  and  0  and  their  derivatives  were  set  to  zero  for  rj  >  o  except  the 
values  for  f  were  taken  to  be  the  value  at  the  old  ^  extended  out  to  the  new  t]c>a. 


The  adjustment  for  mixed  convection  is  slightly  more  complicated.  Again  all  values  were 

set  to  zero  except  the  values  for  f,  f ',  (-^-)0,  (-^-)0,  f0  and  PD.  They  were  adjusted  as 

V  H  }  \  d£  ) 

follows. 


f (J0C,  =  (f,  -  fo 

(lr)°' = “2z°n°(n + IK"’  = fo<  +  t  p°'  +  p- 


■o  i)l*Z 

q 


(45) 


This  adjustment  introduced  less  than  1.0  %  difference  in  the  local  Sussclt  number  for  I’r-  0.7 
as  compared  to  a  constant  run  at  rj  =  30.  'ITtis  error  was  due  primarily  to  r being  less  than  30 
initially  rather  than  to  the  adjustment  itself. 


It  was  found  that  as  <*;  was  increased,  errors  resulted  owing  to  the  growing  boundary  layer 
thicknesses.  By  using  Ayj  =  0.01  numerical  errors  were  reduced  for  calculations  at  high  values  of 


jv:  vj  v.'  wTftw.Tr.'irev'w’i^Yi 


the  curvature  parameter  £.  It  was  found  that  the  step  size  Af  did  not  affect  the  solution 
appreciably  and  A£  =  0. 1  was  used. 


RESULTS  AND  DISCUSSION 


Numerical  results  were  obtained  for  the  case  of  buoyancy  assisting  flow.  They  cover 
Prandtl  numbers  of  0.1,  0.7,  7,  and  100;  values  of  D0  of  0,  0.02,  0.1,  0.5,  1,  2  as  well  as  the  pure 
free  convection  case;  for  a  power  law  temperature  distribution  of  the  form  Tw  —  T„  =  ax",  with 
the  exponent  n  varying  in  the  range  from  -0.4  <  n  <  0.5  for  mixed  convection  and 
— 0.5<n<0.5  for  pure  free  convection.  These  ranges  are  within  the  physical  limits  of  n  as 
determined  in  a  manner  outlined  by  Gebhart  [21].  For  pure  forced  convection  the  limits  are 
—0.5  <  n  and  for  pure  free  convection  the  physical  limits  are  -0.6  <  n  <  1.0.  The  case  of  £20  =  0 
corresponds  to  pure  forced  convection  and  Q0  =  oo  corresponds  to  either  uixs  =  0  (pure  free 
convection)  or  rQ  =  oo  (the  flat  plate  solution  in  mixed  convection).  The  solution  method 
involves  a  system  of  simultaneous  equations  and  the  results  depend  heavily  upon  the  choice  of 
rj„.  References  [4]  and  [13]  incompletely  attributed  the  errors  at  high  curvatures  to  a  large  step 
size  in  the  radial  direction.  It  was  found  that  decreasing  the  step  size  does  cause  a  corresponding 
decrease  and  a  better  accuracy  in  the  local  Nusselt  number  at  high  curvatures,  but  that  the 
velocity  profile,  for  pure  free  convection  in  particular,  was  more  dependent  on  the  choice  of  >7^. 

First,  the  results  for  pure  free  convection  will  be  presented.  The  results  of  Lee  et  al.  [13] 
for  the  UWT  case  were  recalculated  using  a  smaller  step  size  of  Arj*  =  0.01  and  a  variable 
which  was  allowed  to  increase  to  45.  The  radial  coordinate  rj*  used  in  [13]  is  related  to  >7  by 
17*  =  >;/s/2~.  As  an  example,  for  Pr  =  0.7,  present  calculations  for  AN=  50  using  Arj*^  =  0.035 
and  »?*„„  =  7.071  as  was  done  in  [4]  yield  the  result  Ni^Gr,,^4  =  16.1,  whereas  calculations  using 
Ay*  =  0.0125  and  =  30,  as  was  done  in  [13],  gives  Nu,Gt;i/4  =  6.93.  For  Arj*  =  0.01  and 
►j400  =  45,  the  present  calculation  gives  Nt^Grj114  =  6.46.  This  is  only  7.3  %  lower  than  the 
result  of  [13],  and  further  reductions  in  Ay*  or  increases  in  did  not  appreciably  change  the 
local  Nusselt  number.  This  led  to  the  choice  of  Arj*  =  0.01  and  a  variable  up  to  45. 
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llie  Nu,(}r,  l/4  results  are  listed  in  Table  1.  It  demonstrates  that  the  local  Nussclt  number 
increases  with  increasing  curvature  and  increasing  Prandtl  number  for  all  n  values  investigated. 
It  also  shows  that  as  the  curvature  increases  the  effect  of  Pr  decreases  and  the  value  of  the  local 
Nusselt  number  converges  to  an  asymptotic  value.  Table  2  shows  the  Nig/Nu,  uwx  values  as  a 
function  of  the  curvature.  These  values  approach  one  as  the  curvature  increases,  independent  of 
n  for  all  Prandtl  numbers.  This  is  anticipated  because  as  the  curvature  increases  the  first  two 
terms  in  the  energy  equation  become  dominant  and  the  equation  becomes  independent  of  both 
Pr  and  n.  The  average  Nusselt  number  results,  NuL(ir, 1 4,  are  shown  in  fable  3.  They  show 
trends  similar  to  that  for  the  local  Nusselt  number. 

'lire  effect  of  surface  curvature  on  the  local  Nusselt  number  is  illustrated  in  I  ig.  1  for  the 
IJWT  case.  It  is  seen  that  as  the  curvature  increases,  the  local  Nusselt  number  increases  at  a 
faster  relative  rate  for  lower  Prandtl  numbers,  figure  2  illustrates  the  effect  of  the  exponent  n 
on  the  local  Nusselt  number  for  the  flat  plate.  It  is  seen  that  the  local  Nusselt  number  increases 
with  increasing  n.  It  also  shows  a  slight  dependence  on  the  Prandtl  number  such  that  the  higher 
the  Prandtl  number  the  less  is  the  deviation  from  the  UWT  case. 

Results  for  P'(f ,  0)  are  given  in  Table  4.  It  is  noted  that  the  general  trend  is  that  P'(£,  0) 
increases  with  increasing  curvature,  as  noted  by  Minkowycz  and  Sparrow  [12].  Exceptions  to 
this  may  be  found  at  high  curvatures  and  may  be  due  to  numerical  errors.  It  decreases  with 
increasing  n  and  increasing  Prandtl  number  as  noted  by  Chen  et  al.  [22]  for  the  vertical  plate. 
'Hie  relationshop  between  P'({,0)  and  n  is  opposite  to  that  for  the  local  Nusselt  number.  This 
can  be  explained  by  realizing  that  lower  Prandtl  numbers  will  result  in  larger  velocity  gradients 
at  the  wall  whereas  larger  Prandtl  numbers  will  yield  a  larger  wall  temperature  gradient  and 
therefore  a  larger  Nusselt  number. 

Next,  the  results  for  pure  forced  convection  will  be  presented.  Results  of  I ,ee  et  al.  [4]  for 
the  UWT  case  were  found  to  be  in  error  at  high  curvatures  for  the  same  reasons  as  mentioned 
in  the  case  of  pure  free  convection.  As  an  example,  results  in  [4]  were  obtained  using 
A»j  =  0.05  and  =  10.  At  AF  =  50  for  Pr  =  0.7  and  n  =  0,  this  gives  Nu,Rex1/2  =  16.2,  which 
compares  with  Nu,Re,-|/J  =  6.72  for  Ar;  =  0.01  and  -  45.  further  decreases  in  Ary  and  rj „  did 
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not  appreciably  lower  the  value  of  NuxRex  1 2,  so  the  problem  was  run  using  Arj  =  0.01  and  ^ 
was  varied  up  to  45. 

With  increasing  for  increasing  curvature,  it  is  of  interest  to  determine  the  degree  to 
which  the  boundary  layer  assumptions  remain  valid.  For  the  boundary  layer  theory  to  be  valid, 
(r4  —  r0)/x  must  be  of  the  order  of  magnitude  of  the  term  (ReJ/2  +  GrJ'4) -l  where  (r4  —  rQ)  is  the 
thickness  of  the  boundary  layer.  The  equations  for  A  and  n  tan  be  combined  to  give 


-(Re 


1/2  4  <«/')=  -*-[(!  +  Ar,d)'12-  !] 


(46) 


Typically,  »;4  is  defined  as  the  value  of  »j  where  f  and  4>  are  within  0. 1  %  of  their  maximum 
value  inside  the  boundary  layer.  For  A  =  1,  rjs  is  about  20  and  the  right-hand  side  of  equation 
(46)  is  about  7.  If  [(r4  —  r0)/x](ReJ/2  +  Gr]/4)  remains  constant  and  A  increases  to  50,  rj6  will 
increase  to  about  620.  These  values  show  that  the  range  of  can  extend  very  large  as  the 
curvature  increases  and  the  boundary  layer  assumptions  will  remain  valid.  It  also  verifies  the 
analysis  presented  in  [13]  and  [12]  that  for  free  convection  under  the  UWT  condition,  the  term 
that  is  a  measure  of  the  boundary  layer  thickness  may  actually  decrease  despite  an  increase  in 
ij4.  This  would  happen  in  the  example  if  rjs  was  less  than  620  but  greater  than  20. 

The  Nu,Re, 1/2  results  for  pure  forced  convection  are  listed  in  fable  5.  The  table  shows 
that  as  the  curvature  increases  the  local  Nusselt  number  increases  lor  all  Pr  and  n  that  were 
investigated.  It  also  shows  that  as  the  curvature  increases,  the  effect  of  Pr  diminishes  and  the 
Nusselt  number  converges  to  an  asymptotic  value.  Included  in  Table  5  are  the  f  '(£ ,  0)  results. 
It  should  be  noted  that  for  pure  forced  convection  the  momentum  equation  is  not  coupled  to 
the  energy  equation  and  thus  P'(£,  0)  is  independent  of  n  and  Pr.  It  is  seen  that  as  the  curvature 
increases  P'(£,  0)  increases,  as  has  been  noted  in  other  studies  (for  example,  Seban  and  Bond 
[23]).  In  Table  6  are  listed  the  NuJNu^,,^  results  as  a  function  of  n.  The  effect  of  n  is  seen  to 
diminish  as  the  curvature  increases.  Results  for  the  average  Nusselt  number  follow  the  same 
trends  as  the  local  Nusselt  number,  as  can  be  seen  in  fable  7. 
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Figure  3  is  a  plot  of  NuJNu^  vs.  AF  and  shows  that  the  local  Nusselt  number  ratio 
increases  faster  with  increasing  curvature  for  lower  Prandtl  numbers  and  becomes  nearly  linear 
at  higher  curvatures.  This  trend  is  similar  to  that  for  natural  convection.  The  effect  of  the 
exponent  n  on  the  flat  plate  is  illustrated  in  Fig.  4.  It  shows  that  as  n  increases  the  local  Nusselt 
number  increases  and  that  as  I’r  increases  the  deviation  from  the  UWT  becomes  less. 

Calculations  for  mixed  convection  were  done  using  Ay  =  0.01  and  was  varied  up  to  45. 
The  independent  variable  £  was  varied  from  0  (corresponding  to  pure  forced  convection  for  the 
vertical  flat  plate,  y  =  1  and  A  =  0)  to  the  point  where  either  A  >  50  or  y<  0.1.  Tables  8 
through  13  present  the  results  for  Nu*/(ReJ/J  +  GrJ/4)  and  Tables  14  through  19  present  the 
results  for  0).  Tables  8  and  14  present  the  special  case  of  mixed  convection  for  vertical 
plates  (A  =  0,  Q0  =  oo).  For  mixed  convection  along  vertical  cylinders,  although  calculations 
were  done  using  the  parameter  ila  to  combine  the  curvature  and  buoyancy  parameters  in  one 
independent  variable  £,  results  are  more  meaningful  when  presented  in  terms  of  the  curvature 
parameter  A  and  the  mixed  convection  parameter  y.  This  is  shown  for  the  UWT  case  in  Figs.  5 
through  8  for  the  local  Nusselt  number.  In  these  figures  the  solid  lines  represent  the  actual  runs 
for  given  values  of  while  the  dashed  lines  are  found  by  interpolation  of  the  computer  output 
for  the  values  of  A  presented.  In  all  cases  the  curve  for  the  flat  plate  (A  =  0)  is  concave  upward 
and  as  x  decreases  from  1  to  0  along  the  constant  curvature  curve  of  A  =  0,  the  value  of 
Nu,/(ReJ/J  +  GrJ/4)  decreases  to  a  minimum  in  the  vicinity  of  y  =  0.6  then  increases  to  its  pure 
free  convection  value  at  y  =  0.  This  is  seen  more  clearly  in  Fig.  9  for  the  flat  plate, which  also 
demonstrates  that  for  n  #  0,  the  curves  run  essentially  parallel  to  the  UWT  curve.  Figures  5 
through  8  clearly  show  that  as  the  curvature  increases  for  constant  values  of  y  the  local  Nusselt 
number  increases.  The  results  shown  in  Tables  8  through  13  thus  reflect  the  combined  effects  of 
A  and  y  on  the  Nusselt  number.  That  is,  if  the  Nu,  value  initially  decreases  with  increasing  £, 
as  it  does  for  the  higher  values  of  Q„,  then  the  effects  of  mixed  convection  are  dominant.  If  it 
initially  increases  with  increasing  £,  then  the  effects  of  curvature  are  dominant. 

To  observe  the  effect  of  O0  on  the  Nusselt  number  results,  one  refers  to  Figs.  5  through  8 
and  notices  that  when  £10  -*  oo  cither  the  flat  plate  solution  is  approached  or  the  pure  free 


convection  case  is  approached.  If  the  curves  cut  off  at  /  <  0. 1  were  extended  they  would  rise 
sharply  to  run  between  the  curves  for  Q0  =  0.5  and  ft0  =  oo,  the  pure  free  convection  case. 

The  local  Nusselt  number  results  obtained  at  high  curvatures  in  this  study  were  much 
lower  than  those  obtained  by  Lee  et  al.  [4],  as  has  already  been  pointed  out.  It  was  found  that 
results  for  higher  Prandtl  numbers  were  affected  less  by  a  decrease  in  the  step  si/e  and  an 
increase  This  perhaps  may  be  due  to  an  insufficient  value  of  >j^,  but  it  is  most  likely  due  to 
the  fact  that  for  lower  Prandtl  numbers,  the  first  two  terms  of  the  energy  equation  essentially 
dominate  the  equation,  whereas  for  higher  I’randtl  numbers  other  terms  in  the  energy  equation 
remain  significant.  I,ee  et  al.  [4]  made  the  statement  that  for  A  >4  the  solution  for  mixed 
convection  was  independent  of  y,  Pr  and  n.  Although  this  is  true  as  A  -»  oo,  and  appears  nearly 
true  for  the  results  of  [4]  it  is  not  the  case  for  the  present  results.  Lor  A  <  50,  the  effect  of 
varying  Pr  from  0.1  to  100  remains  significant  and  the  pure  free  convection  result  Nu^Gr,1'4  does 
not  equal  the  pure  forced  convection  result  Nig.Re,,1  2  which  fee's  statement  seems  to  imply. 
For  low  Prandtl  numbers  this  is  nearly  the  case  but  for  higher  Prandtl  numbers  this  is  not  the 
case. 

The  average  Nusselt  number  results  were  also  obtained.  Those  for  the  UWT  case  arc- 
presented  in  Figs.  10  through  13  and  those  for  the  flat  plate  with  nonuniform  T„(x)  are 
illustrated  in  Fig.  14  ITie  trends  of  the  curves  are  similar  to  those  for  the  local  Nusselt  number 
The  value  of  P'(£,0),  listed  in  Fables  14  through  19,  in  general  decreases  initially  and  then 
increases  as  £  increases.  This  is  similar  to  the  results  of  the  Nusselt  number  but  the  magnitude 
of  the  increase  with  increasing  curvature  is  less  pronounced  for  P'(£,  0). 

For  practical  purposes,  correlation  equations  were  developed  for  the  local  and  average 
Nusselt  numbers.  For  pure  free  convection,  for  0.1<Pr<100,  —  0.5  <  n  <  0.5,  and 
0  <  A  <  50  (A  =  An),  the  correlation  equation  for  the  local  Nusselt  number  is  given  by 

NuxGr;,/4  =  *N(Pr)[AN(A)  +  f,.N(Pr)A](l  f  VNWN)  (47) 


where 


<xN(Pr)  =  0.75Pr  '  [2.5(1  +  2Pr'‘  +  2Pr)] 


(1-0  7-') 


An(A)=-  1.4' 


f,  N(Pr)  =  0.045  i  0.27Pr~ 


VN  =  {[0.88  f  0.5  cxp(  -  2Pr  ■  )]  -  0.7^n}n 


WN  =  exp{  -  [0.07  +  7.5exp(  -  3.6Prl/IO)]A°  7} 


The  average  Nussclt  numbers  arc  correlated  bv 


NuLGr,.1/4  =  i-aN(I»r)[BN(A)  +  f2  N(Pr)A](l  +  VN) 


where 


(1-  0.7SA) 


Bn(A)=  1.4' 


f2  N(Pr)  =  0.03  +  0.21  Pr“ 


V  =  (3VnWn  — nW)/(3  +  nW) 


W  -  exp(  -0.1  IA  2) 


In  equations  (53)  -  (57),  A  now  stands  for  AN  with  x=  L. 


For  the  case  of  pure  forced  convection,  for  0. l<Pr<)00,  — 0.4<n<0.5,  and 
0  <  A  <  50  (A  =  Af),  the  correlation  equations  are 

NuxRc-,/2  =  aK(Pr)[AK(A)  +  f,,F(Pr)A](l  +  VFWr.)  (58) 


where 


■e.  -V.  -  o  '  W  V  o  -  '  «.■  hJ-  ■ 


(59) 


otF(Pr)  =  0.339Pr,/3[  I  +  (0.0468/ Pr)2/3j“ 1 /4 


A,.(A)  -  (I  f  0.4AI/2) 


m 


f,  ,.(Pr)  =  (0  04  +  ().3Pr  0  J5) 


(61) 


VF-  =  {[1.17  +  1  !cxp(  -5  6Prl/10)J  -  0.92n}n 


(62) 


WF  =  exp[  -  (0.07  f  0.3Pr",)  (2)A3/5] 


(63) 


and 

Nu,  Rcl1/2  =  2xp(Pr)[Ap(A)  +  f2,,.(Pr)A](l  +  Vp)  (64) 

where 

Bp(A)  =  (1  +  0.25AI/2)  (65) 

f2  p(Pr)  =  (0.025  -P  0. 1 5Pr-0  45)  (66) 

VF  =  Vp-expf  -  (0.06  +  0.17Prl/3)A0  65]  (67) 

In  equations  (64)  -  (67),  A  stands  for  AF  with  x=  L. 

It  should  be  noted  that  the  form  of  these  correlation  equations  is  such  that  for  the  UWT 
case  VN,  VK,  VN,  Vh  are  zero  and  the  last  terms  in  equations  (47),  (53),  (58),  and  (64)  can  be 
omitted.  Also  the  form  of  AN,  Ah,  BN  and  Bh  are  such  that  for  A  =  0  they  become  one  and  for 
the  flat  plate  UWT  case  only  aN  and  ay  are  needed.  The  expression  for  aN  is  taken  from  I  de 
[24]  and  the  expression  for  xh  is  taken  from  Churchill  and  Ozoc  [25],  'Hie  error  in  equations 
(47)  and  (53)  is  less  than  8.5  %  for  the  UWT  case  and  less  than  13.5  %  for  -  0.5  <  n  <  0  5. 
The  maximum  error  in  equations  (58)  and  (64)  is  about  9.2  %  for  the  UWT  case  and  about 


Following  Churchill  [20],  the  correlation  equation  for  Nusselt  numbers  in  mixed 


convection  is  expressed  by  the  form 


This  form  of  correlation  has  been  found  to  give  an  accuracy  of  about  5%  for  0. 1  <  Pr  <  100  for 
flat  plates  with  m  =  3  and  was  also  verified  experimentally  for  air  in  the  UWT  case  [14,  15]. 
For  the  present  study  with  a  single  mixed  convection  parameter  x,  the  corresponding  correlation 
equation  can  be  represented  by 


This  equation  was  found  to  be  valid  for  the  case  of  flow  along  a  vertical  cylinder  with  a  power 
law  variation  in  the  surface  temperature  in  the  following  manner.  It  was  found  that  as  A 
increases  the  corresponding  value  for  m  decreases  to  1.  'I'his  can  be  seen  by  observing  the 
decrease  in  the  concavity  of  the  dashed  lines  in  the  Nusselt  number  figures  as  A  increases.  It  is 
also  seen  that  the  concavity  disappears  more  quickly  for  lower  Prandtl  numbers.  The  graphical 
technique  of  Churchill  [20]  was  used  to  determine  the  relationship  among  m,  A,  and  Pr  for  the 
UWT  case.  Figure  15  shows  a  sample  graph  for  Pr=0.7.  In  this  graph  NuF  is  the  Nusselt 
number  as  a  function  of  AK  =  A,  Pr,  and  n  for  the  pure  forced  convection  case  and  Nun  is  the 
Nusselt  number  as  a  function  of  AN  =  A,  Pr,  and  n  for  the  pure  free  convection  case.  That  is, 
they  are  calculated  from  the  endpoints  of  the  curves  in  the  Nusselt  number  figures  at 
X  =  0  and  1  for  lines  of  constant  A  for  given  Pr  and  n.  Reference  lines  for  various  values  of  m 
are  plotted  and  points  corresponding  to  various  curvatures  are  marked  with  different  symbols. 
For  Pr  =  0.7  it  is  seen  that  as  the  curvature  increases  from  A  =  0  to  50  the  value  of  in  decreases 
from  about  3  to  about  1.  T  his  was  done  for  the  other  Prandtl  numbers  lor  the  UWT  case  for 
both  the  local  and  average  Nusselt  numbers  and  the  relationship 


was  developed.  It  was  found  that  equation  (69)  was  also  valid  for  the  flat  plate  for 
-0.4<n<0.5  with  m  =  3  and  the  values  of  Nu^Re;1'2  and  Nu/jt;"4  were  those  obtained  for 
the  given  values  of  Pr  and  n,  with  A  =  0.  For  A  >  0,  it  was  found  that  equation  (70)  could  be 
used  for  n  >  0  but  was  inconsistent  for  n  <  0.  I  lowever,  a  better  relation  for  n  <  0  could  not  be 
found.  For  the  ranges  of  Pr,  A,  and  n  >  0  studied,  equation  (69)  along  with  the  m  expression 
from  equation  (70)  produce  results  valid  to  within  about  9%  when  using  computer  output  for 
Nuf  and  Nun  and  to  within  about  14%  when  using  equations  (47)  and  (58).  The 
corresponding  correlation  equation  for  the  average  Nusselt  number,  where  Nu„,  Re„,  Gr„,  and  x 
are  replaced,  respectively,  by  N'uL,  ReL  GrL  and  xL,  was  also  found  to  be  valid  with  the  same 
accuracy  for  n  >  0.  A  reliable  relationship  for  n  <  0  was  not  obtained,  but  equations  (69)  and 
(70)  usually  will  provide  a  satisfactory  result  for  both  the  local  and  average  Nusselt  numbers. 

These  correlations  presented  are  in  a  different  form  than  those  proposed  by  Lee  et  al.  [4] 
and  Ixe  et  al.  [13].  Their  forms  were  derived  based  on  the  asymptotic  solution  as  A  -* oo.  As 
a  result,  they  proposed  a  logarithmic  correlation  form  for  the  cases  of  pure  forced  and  pure  free 
convection,  which  may  be  difficult  for  use  in  practice.  The  correlations  in  the  present  study  arc 
based  on  the  flat  plate  solution  and  take  advantage  of  the  nearly  linear  relation  between  the 
Nusselt  number  and  curvature  as  shown  in  Fig.  1  and  3 

The  correlation  proposed  in  [4]  for  mixed  convection  is  also  based  on  the  asymptotic 
solution  as  A  -*  oo.  The  correlation  presented  here  is  based  on  the  correlation  for  the  vertical 
flat  plate  proposed  by  Churchill  [20],  which  has  been  verified  experimentally  for  air  [15].  It  also 
is  based  on  the  curvature  parameter  A,  whereas  in  Ref.  [4]  AN  and  AF  were  used  in  the 
calculations  of  the  normalized  Nusselt  number  for  mixed  convection.  As  x  approaches  0  or  1, 
these  calculations  become  unwieldy.  Therefore,  it  is  suggested  that  the  Nusselt  number 
calculations  for  the  pure  free  or  pure  forced  convection,  which  arc  used  in  the  mixed  convection 
calculations,  be  based  on  A. 


CONCLUSIONS 


Mixed  convection  along  a  vertical  cylinder  with  a  power  law  variation  in  the  wall 
temperature  is  investigated  for  the  entire  regime  ranging  from  pure  free  convection  to  pure 
forced  convection.  Results  were  presented  for  the  local  and  average  Nusselt  numbers  and  values 
of  P'(<L  0)  were  tabulated.  Correlation  equations  for  the  local  and  average  Nusselt  numbers  are 
also  presented.  It  is  found  that  the  local  Nusselt  number  NuJfRey2  +  GrJ/4)  increases  with 
increasing  Prandtl  number,  increasing  value  of  the  exponent  n,  and  increasing  curvature  for  the 
entire  mixed  convection  regime  ranging  from  pure  free  convection  (x  —  0)  to  pure  forced 
convection  (x  =  1)-  For  the  vertical  flat  plate,  this  quantity  initially  decreases  as  x  is  decreased 
from  1  and  then  increases  to  the  pure  free  convection  value  as  y  approaches  0.  However,  this 
trend  tends  to  become  a  linear  relationship  as  the  curvature  increases;  that  is,  the  value  of 
Nux/(ReJ/2  +  GrJ/4)  for  constant  curvature  A  will  lie  on  a  line  with  end  points  at  x  =  0  and  x  =  1 
for  that  curvature.  This  trend  is  approached  more  rapidly  for  lower  Prandtl  numbers.  The 
quantity  NuL/(Re{/2  +  Grj/4)  follows  a  similar  pattern  as  that  for  the  local  Nusselt  number.  The 
effects  of  both  the  Prandtl  number  and  n  on  the  Nusselt  number  diminish  as  the  curvature 
increases,  but  remain  significant  for  curvatures  A  <  50  when  Prandtl  numbers  are  large. 
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Variation  of  Surface  Temperature,  Mixed  Convection,  Along  a  Vertical  Flat  Plate 
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IV.  RECOMMENDATIONS  FOR  FURTHER  STUDY 


From  the  present  study  it  was  found  that  there  exists  a  need  to  develop  an  accurate  way  of 
extending  in  the  numerical  computations,  l^crge  values  of  ^  are  required  at  high  curvatures, 
but  starting  out  at  a  large  value  is  both  numerically  inefficient  and  can  result  in  calculations 
beyond  the  precision  capabilities  of  the  computer.  Therefore  tj^  must  be  increased  gradually 
from  its  initial  value  of,  say,  10.  Secondly,  the  decision  must  be  made  as  to  when  *7 „  is 
sufficient.  This  could  be  done  by  comparing  both  the  boundary  conditions  and  some  select 
points  on  the  velocity  and  temperature  profiles  for  the  solution  at  one  value  of  ^  with  the 
solution  for  an  increased  value  of  while  keeping  the  step  si/e  Arj  fixed.  When  there  is  no 
longer  a  change  in  these  values,  then  the  solution  at  that  node  is  accepted  and  the  problem 
advances  to  the  next  node  in  the  £  direction.  Limitations  to  this  might  be  the  memory 
capabiltiy  of  the  computer  or  precision  difficulties  if  the  thermal  boundary  layer  is  vastly 
different  from  the  flow  boundary  layer.  This  can  occur  for  both  low  and  high  Prandtl  numbers. 
The  memory  problem  can  be  helped  by  trying  a  step  size  initially  of  0.05.  With  the  increased 
value  of  17^  this  may  provide  sufficient  accuracy  and  should  reduce  the  number  of  nodes  in  the 
radial  direction. 

Experimental  results  are  needed  for  flow  along  vertical  cylinders  for  the  full  range  of  mixed 
convection,  ranging  from  pure  forced  convection  to  pure  free  convection,  under  various  surface 
heating  conditions.  Experimental  results  are  essential  to  determine  the  extent  to  which  the 
imposed  surface  temperature  variation  in  the  analysis  can  be  realized  in  the  physical  world.  The 
singularity  at  x=0  for  the  case  of  n<0  in  the  power  law  variation  of  the  surface  temperature, 
Tw(x)  —  T„  =  ax",  for  example,  must  be  accounted  for  and  an  experiment  is  the  only  real  way 
this  can  be  checked. 

Variations  of  the  problem  to  include  transition  to  turbulence,  or  unsteady  variations  in 
wall  temperature  could  be  addressed.  These  should  be  done  in  conjunction  with  experiments 
when  possible  to  prove  their  applicability  or  to  provide  information  such  as  when  and  what  type 
of  transitions  to  turbulence  occurs. 
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APPENDIX  A 


TRANSFORMATION  OF  THE  GOVERNING  SYSTEM  OF  EQUATIONS  IN 

CHAP  TER  2 


In  Chapter  2  the  solution  for  the  case  of  natural  convection  along  a  slender,  vertical 
cylinder  with  variable  surface  heat  flux  is  presented.  Consider  a  semi-infinite,  vertical  cylinder 
with  radius  r0  that  is  aligned  in  a  quiescent  ambient  fluid  at  temperature  l'rxj.  The  axial 
coordinate  x  is  measured  upward  when  q„  >  0  and  downward  when  qw  <  0.  The  radial 
coordinate  r  is  measured  from  the  axis  of  the  cylinder.  Ihe  surface  of  the  cylinder  is  subjected 
to  an  arbitrary  heat  flux  q„(x)  and  the  gravitational  acceleration  g  is  acting  downward.  The  fluid 
properties  are  assumed  to  be  constant  except  for  variations  in  density  which  induce  the 
buoyancy  force.  By  employing  the  laminar  boundary  layer  assumptions  and  t.  aking  use  of  the 
Boussinesq  approximation  the  governing  conservation  equations  can  be  written  as: 

Continuity: 

|r(n,|  +  ^n,)-0  M-'> 


Momentum: 


uiu  +vj|t 

ox  or 


-Tfr(r£-)+*,fr-1- 


Energy: 


dx  dr  T  dr  \  dr  / 


(A -2) 


(A  -  .1) 


In  these  equations  u  and  v  are  the  velocity  components  in  the  x  and  r  directions 


respectively;  T  is  the  fluid  temperature;  and  v,  (i  and  a  are,  respectively,  the  kinematic  viscosity, 
the  volumetric  coefficient  of  thermal  expansion  and  the  thermal  diffusivity  of  the  fluid. 
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ITte  boundary  conditions  are 


u  =  v  =  0  , 


-  qwW 


at  r  =  r„ 


u  -*  0  ,  T  —  As  r 


u  =  0  ,  I  =  1^  At  x  =  0,  r  >  rD 


(A -4) 


(A-  5) 


(A  -  6) 


In  writing  the  last  boundary  condition  it  is  assumed  that  the  flow  and  thermal  boundary 
layer  thicknesses  are  zero  at  the  leading  edge  of  the  surface. 

I  he  conservation  equations,  along  with  the  boundary  conditions,  are  transformed  into  a 
dimensionless  form  by  introducing  the  following  dimensionless  variables: 


2  2 

=  ~[/~2 r'°--(Gr**/5)1/s,  * = Y^v/sr''5 


(A-  7) 


f(A,  n)  =  ^(x,r)/[5vr0(Gr//5)1/5],  0(3, ,,)  =  ----  (A  -  8) 

qw(x)  x/k 


=  gfaw  (x)x  /kv2 


(,1  -  9) 


where  »/  is  the  pseudo-similartiy  variable,  ((X,rj)  is  the  reduced  stream  function,  0(X,  tj)  is  the 
dimensionless  temperature,  X  is  the  curvature  parameter,  and  i^(x,  r)  is  the  stream  function  that 
satisfies  the  continuity  equation,  with  u  =  (dik/dr)/r  and  v  =  -  (<?i£/r?x)/r.  The  transformation 


u  =  5^-(Gr//5)2/V(^,  >,) 


(A  -  10) 


I# 


/Vv'N 


8$ 

m* 


gif 

P 

Em 


TO 


■ 


tf!* 


1; 


V.V.V.V.* 
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+ 


f[l/fk,^i^s  +  4,5i-|'sqL's‘ 


V  *  *  ■  / 


dx 


5  kS 


+  f  [2/Sq;J'5-^,3'5  + 


=  5v(GrJC*/5)3/sx_2-p-f"(i,  0 


dr 


(/l  -  12) 


(A  -  13) 


d_ 

dr 


10v(Grx*/5)3/5x-2-r-P'  +  5v(Grx*/5)4/sx 

ro 


4/S.-3  r_p„ 
,2 


(.4  -  14) 


T  -  To,  =  — (GrJ[*/5)_1/50(A,  ,)  (/I  -  15) 

JT  =  j/  gP  Vi/sf/ag  ai_  0,^LVg4^5x|/5') 

dx  *  \  5£v2  /  1\d^  dx  dx  A  w  ' 

V  /  (/I -16) 


3T 

dr 


K  r0 


M  -17) 


(Or// 5) 


'  /5 _T _ pi! 


(/I  -  IS) 


in  which  the  primes  denote  partial  differentiation  with  respect  to  rj. 


Substitution  of  equations  (A- 10)  through  (A- 18)  into  equations  (A - 1 )  through  (A-6) 


results  in 


(1  +  >;A)P"  +  AP'  +  (y  +  4)ff"  -  (2 y  +  3)P2  +  0  =  A(y  -  \)( P'-|^  -  P-^-)  (-4  -  19) 

\  aX  aX  ] 

(1  +  r,X)0"  +  X0'  +  Pr(y  +  4)f0'  -  Pr(4y  +  !)P0  =  PrA(y  -  1  -  P^?-)  (-4  -  20) 

\  aX  ok  J 


({X,  0)  =  P(A,  0)  =  0,  0'(X,  0)  =  -  1 

P(A,  oo)  =  0,  0(A,oo)  =  0 


(-4  —  21) 


where 


flw  dx 


(A  -  22) 


The  system  of  equations  (A- 19)  through  (A-22)  represents  the  general  form  of  the 
transformed  boundary  layer  equations  for  variable  heat  flux,  q„(x),  along  vertical  cylinders  in 
natural  convection.  For  the  case  of  power  law  heat  flux  distribution  qw(x)  =  a\n  ,  one  has  from 


equation  (A-22)  that 


(A  -  23) 


For  long  slender  cylinders,  the  curvature  parameter  X  —  2~(Gtj*/5)  1/5  can  be  large.  To 

*0 

lower  the  maximum  value  of  calculations  in  the  x  or  X  coordinate,  one  introduces  a  new  <J(x) 


variable  defined  by 


A  =  C  i*2 


(A  -  24) 


with  C  =  2 . 5,/s. 


Substituting  equations  (A-23)  and  (A-24)  into  equations  (A- 19)  through  (A-22)  results  in: 


Momentum: 


(1  +  a ,r,)P"  +  a,P'  +  a2fP'  +  a3P2  4-  a 40  =  a/p'-|^  -  P-f^ 


(A  -  25) 


(1  4-  a,ij )0”  4-  a ,0'  4-  Pr  a2f0'  4-  Pr  a6P0  =  Pr  a5f  . 0'-^- 


Boundary  Conditions  : 


f(i,0)  =  P(£,0)  =  o,  £?'(<?.  o)  =  -  l 

Ptf .  oo)  =  0(f.  oo)  =  0 


(26) 


(27) 


where 

a,  =  Ci*2  ,  a2  =  n  4  4  ,  a?  =  -  (2n  4-  3) 
a4  =  1  ,  a5  =  {(n  -  l)/2  ,  a6  =  -  (4n  +  1) 


(28) 


The  proposed  method  of  solution  ,  developed  by  Ixe  et  al.  (1986b)  requires  boundary 
conditions  in  terms  of  f,  f  and  0  at  ij  =  0.  In  equation  (27)  one  has  Q'(/.,  0)  =  -  1.  This 
necessitates  an  additional  transformation.  Let 


*«,»»)  =  <?«,. i)/0((,0) 

(29) 

Then: 

=  0) 

(30) 

(31) 

(32) 

Substituting  0'({,  0)  = 

-  1  from  equation  (27)  into  equation  (31)  for  r\  =  0  results  in 

0(1,  0)  -  “  1 

(33) 

Substituting  equation  (33)  into  equations  (30)  through  (32)  leads  to 


Olt.n)  =  -*«,.»)/*'«.  0) 


(34) 


■fc'A'.V.fi 


l  |.l  Kkkf  It'.li'li*  tJ  «.l  «.<  M  ».»  »*« 


0'(£,  »?)  =  —  4>’(Z,  *i)l<t>'{i<  0) 


a  l~  4>(£,>?)  ~ 
^  L^«.o)_ 


35) 


M  -  36) 


(/l  -  37) 


Substituting  equations  (A-29)  through  (A-37)  into  equations  (A-25)  through  (A-28)  results  in  : 

(1  +  a ,,)P"  4-  a,r  +  a2fP'  +  a/2  +  a7<*>  =  a5(f^|-  +  P-^- j  (/<  -  38) 


(1  +  +  a,^'  +  a2Pr  f<t>'  +  a6Pr  P<£ 


=  Pra5  f-^-r  f(U)  Ir- 


(A  -  39) 


Boundary  Conditions  : 


f(£,0)  =  P(£,0)  =  0  W,  0)  =  1 


f(£,  oo)  =  <£(<?>  cx>)  =  0 


(A  40) 


where  a,  through  a6  are  as  before  and 


*'«,  0) 


(/I  —  41) 


The  boundary  conditions  arc  now  in  the  proper  form  to  be  used  in  the  finite  difference  solution 
method  described  in  Appendix  B. 

The  physical  quantities  of  interest  are  the  local  Nussclt  number  Nux,  the  average  Nusselt 
number  NuL,  the  local  wall  shear  stress  tw,  the  axial  velocity  distribution  u,  and  the  temperature 
profile  4>(Z,  rj).  The  local  Nusselt  number  is  defmed  by 


A  -  42) 


«v-y'>'AA-V.-v/v'/v:v' 


.*  -* 


It*  i. 


m, 

r 

f 

f 

i* 


m 

m 


m 


m 

r 

V. 


is 
$$$ 


■vrr 

*AW 
.  ^VV- 

-V\- 

»  V* ■  > 

>S;%V 


v%.-  v 

.  *  *  v  * 


where  h  =  qw/(T„  -  T„).  Substituting  qw  from  equation  (A-4)  into  equation  (A-42),  followed 


by  the  transformation,  results  in 


=  — L 


Nux(Cir//5)-"  = 


0({,0) 


(A  -  43) 


Substituting  equation  (A-33)  into  equation  (A-43)  one  obtains 


Nux(Gr//5f1/S  =  -<*>'(£, 0) 


(A  -44) 


The  average  Nussclt  number  is  defined  by  Nu,  =  h~  where  the  average  heat  transfer 


coefficient,  h,  is  obtained  from 


(A  -  45) 


Substituting  h  into  the  definition  of  h  and  carrying  out  the  integration  one  obtains: 

NuL(GrL*/5r,/5  =  -  L~(4+n)/5  fLjc(4+n)/s^'(^  0)x~’dx  (A  -  46) 

•>0 


From  equations  (A-7)  and  (A-24)  one  can  write 


I",  ,  > l  — n  ~ I 
c_  (x/r0)  ,/io 

^  Gr  * 


x  =  r0[Gro^I0]'«'-n) 


x-,dx  =  -ri5_r'd{ 

1  —  n 


(A  -  47) 


(A  -  48) 


(A  -  49) 


where  Gr0*  =  Gr„*  at  x  =  r0;  that  is,  Gr0*  =  gfi(a  r„n  )  r  *lkv2.  Substituting  equations  (A-48)  and 


(A-49)  into  equation  (A-46),  one  gets 


NuL(GrL*/5)'1/5  =  -  fSVtf.  0)di 

1  n  Jn 


(A  -  50) 


4 


where 


r’I.J 


v  vv  V^«  V  V  v.v,^ 

s.y 


<J,  =  ^  at  X  =  I,, 

M  =  (7  +  3n)/(  1  -  n),  N  =  (8  4-  2n)/(n  I ) 


(A  51) 


As  £  -*  0,  equation  (A- 50)  approaches 


Nu,(Gr|*/5) 


•/St"'/* 


4  f n 


— — </>'(().  0) 


(A  52) 


Next,  from  the  definition  of  local  wall  shear  stress 


T«  =  /Jrr 

<;r  I  r=  r„ 


one  obtains,  by  substituting  equation  (A- 13)  into  equation  (A-53), 


(A  S3) 


MS* 

m 

m 


Tw  —  5  — ^<irx*/5)1/sf"(£,  0) 


from  equation  (A- 10),  the  axial  velocity  distribution  can  be  written  as 


(A  S4) 


~r  =  S(Gr//5)2/sf({, ,,) 


and  the  temperature  profile  »j)  is  given  by  equation  (A- 20). 


(I  SS) 
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APPENDIX  B 


METHOD  OF  SOLUTION  FOR  THE  SYSTEM  OF  EQUATIONS  IN  CHAPTER  2 


Equations  (A-38)  through  (A-40)  in  Appendix  A  constitute  a  system  of  nonlinear  partial 
differential  equations  in  the  ({,  tj)  coordinates  with  parameters  Pr  and  n.  The  method  of  I  ee  et 
al.  (1986b)  is  employed  to  solve  this  system  of  equations  for  natural  convection  along  a  vertical 
cylinder  with  q„  =  ax”. 


The  first  step  is  to  convert  the  terms  involving  d/d£  in  the  following  manner. 


—  =  —  (H-H0)-q(— -)0  =  —  (II -Ho) 

d(  A£  v  M  di  A£  0 


(B-  1) 


where  II  is  a  function  of  ({,  rj),  H0  =  1 1„  +  (q/p)A.J;(<?I !/<?{)„,  and  the  subscript  "o'  denotes 
quantities  at  {  -  A£.  The  values  for  p  and  q  arc  p=  1  and  q  =  0  at  {  =  0  and  f  =  A£, 


and  p  =  2  and  q  =  1  thereafter  for  {  >  2A<f.  Thus  we  have 


dt,  M K  o) 


P  /  **({. 0)o  ,  I,  ,  P 


**K.O) 


~  *  )4>  +  ~  ^o) 


(R-2) 


{B-  3) 


(B-  4) 


f0  =  f0  +  -^A{ 


(f> 


fWo  +  ^-A{ 


(«• 


(B~  5) 


{B-  6) 


!;! 

t'tC> 


:S| 

SSl 
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■Sss 
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*• .  •»  • 

•  /-  «  '  *  ■ 

>>> 


•:.v,v 


krta" 


,va; 
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w 
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\-v 


N  -%A*  \ 

. » .  Ai 


■.-■v-iss 


KX^.k  PJ 
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4>'(Z<  0)o  =  4>'(t,  0)  +  ^rAt 


q  A  J  dP(t,  0) 


dt 


(B  -  ?) 


(/?-  8) 


The  f  and  4>  equations  arc  then  linearized  by  using  the  simple  technique  for  the  product  of 
two  arbitrary  functions  F  and  G  ,  defined  as 


FG  =  FG  4-  GF  -  FG 


(B  ~  9) 


where  F  and  G  indicate  the  values  to  be  guessed.  Applying  equation  (11-9)  to  the  terms 
fP' ,  P2 ,  f<f>' ,  and  f <f>  ,  one  obtains 


fP'  =  fP'  4-P'f-P'f 


(//  10) 


P2  =  2PP  -  f 2 


<//-  11) 


ty  =  f<y  +  4>'t-f4>' 


(B-  12) 


V<i>  =  Pqi  +  P<£  -  ftf, 


{B-  1.1) 


Substituting  equations  (B-2)  through  (B-13)  into  equations  (A-38)  through  (A-40)  results 


in: 


Momentum  equation: 

(1  4-  a„)P"  +  [a,  +  (a2  -  a5  p/A{)  f  +  as  P/A<?  fc]  P' 

-I-  [2(a,  4-  a5  p/A£)  P  -  a5  p/A£  PJ  P  4-  [(a2  -  as  p/A£)P']  f  +  a7< />  (/i  -  14) 


=  [(a2  -  a 5  p/A£)fP'  4  (a,  t  a5  p/A^f'2  ] 


Fnergy  equation: 


■2S£vj 

'  I 


W- 


v» 

fcSSq 


4M 


V*V\0 


'Z'/.V, 


£f 


■  fc** 


■fc  *  *  ■ « 


•;:^r 

„  *  ■>  fc,  ■* 


Ilic  system  of  equations  has  now  been  reduced  to  a  system  of  quasi-lincar  ordinary 
differential  equations.  By  defining  g  =  f  and  applying  the  weighting  factors 


[see  Lee  et  al.  (1986b)]  one  can  write  equations  (18)  and  (19)  in  the  finite-difference  form  as: 

f,  -  f,_  i  -  +  g,_ , )  =  0  ( 


Ao/(Aij)  +  a0 g,  +  ot+]g(+1)  +  A2g/  +  A3f,  +  /\44>i  =  A5 


Bq/(Ai;)  +  /?o <t>,  +  P+\4>i+\)  +  B2<0(  +  Bjg,  +  Baf(  =  B5 


where  Ar;  is  the  step  size  in  the  rj  direction  and  the  subscript  "i"  refers  to  values  at  the  nodal 
point  r/,. 

In  equations  (24)  and  (25)  the  weighting  factors  arc  defined  as: 

ot-|  =  VV,(  -  Z(_,/2)  ,  a+1  =  W(<Zi+I/.2)  ,  a0  =  -  a_,-a+| 

0-i  =  -  z*/-i/2) ,  0+i  =  wdz*/+i/2)  -  0o  0-i  _0+i 

(26) 

Wffz)  =  z/(  1  -  exp(  -  z)) 
z  =  Aijfa,  +  A,)/A0  ,  z*  =  Aij(a,  +  B,)/B0 

Equations  (23)  through  (25)  are  applied  to  the  interior  points  On  the  boundaries. 


equation  (16)  is  written  as 


f|=gi=0.  <0i=l 

At] 

—  ^n—  1  3"  8n-l)  —  8p  “  <0n  —  ^ 


where  the  subscript  n  is  the  number  of  nodes  in  the  »j  direction. 


Equations  (23)  through  (27)  constitute  a  system  of  algebraic  equations  that  can  be  written 


in  the  matrix  form 


[A][X]  =  [B] 


where  [A]  is  a  band  matrix  of  order  3n  and  bandwidth  seven.  The  array  [X]  which  contains 
the  solution  in  the  form  (f, ,  g, ,  4>u  f2,  g3,  <j>2,  ....  f„,  g„,  <0„)T  is  a  column  matrix  of  order  3n 


.  *  ■  '.  »  *'*  * 

'  k  S  -  K  • 

~A~ 


/,  v\ 

•-'/-‘.vC* 


m 

m 

■K  -  k  ' 1  C- 

.VMS', 


a  »  »  ■_#»  »  m 

v;,* 

.-.-'.V 

■  «  •  -  ■  a 

sV.s1. 

V^V 

:  s 

•  » « • , 

V  V  V 

y;v;y 

'A/  V 

*VV> 

>vv 


sv*--: 

^%-v: 

£;-Xy 


■■  v:  ,‘.v  '  zXMWi*' 


>>> 

r. 

*\V-\ 


where  [A]  is  a  band  matrix  of  order  3n  and  bandwidth  seven.  The  array  [X]  which  contains 
the  solution  in  the  form  (f, ,  g, ,  f2,  g2,  <f>2,  ...,  fn,  gn,  4>y  is  a  column  matrix  of  order  3n. 

The  matrix  [B]  is  a  column  matrix  of  order  3n  which  contains  the  right  hand  sides  of  equations 
(B-23)  through  (B-25)  and  (B-27).  I'he  matrix  [A]  is  approximately  diagonally  dominant  and 
equation  (B-28)  can  be  solved  by  the  Gaussian  elimination  technique  with  high  accuracy. 

To  obtain  values  for  f"  and  <j>',  the  results  for  P  and  <f>  along  with  the  boundary 
conditions 

P"«,0)  =  -air'«.0)-44  (ft  -  29) 

«T(£.0)=-  a,*'«.0)  (It  -  30) 

P"(£,  oo)  =  4>"(L  oo)  =  0  (It  -  31) 

are  used  in  a  cubic  spline  interpolation  routine  [sec,  for  example,  Burden  and  Faires  (1085)]. 
As  the  solution  converges,  the  boundary  conditions  become  exact  and  the  values  for 
P'  and  <£'  can  be  obtained  with  high  accuracy. 

The  present  method  of  solution  employs  a  quasi-lineari/ation  of  the  original  nonlinear 
system  of  equations  and  requires  initial  guesses  for  f,  P,  P',  <f>,  and  <t>' .  l  he  Hat  plate  solution 
(i.e.,  £  =  0)  for  the  uniform  surface  heat  flux  case  (IJIU;)  for  Pr  =  0.7  was  obtained  and  these 
results  were  used  as  the  initial  guesses  for  all  other  combinations  of  Pr  and  n  at  £  -  0.  For 
£  =  A£,  initial  guesses  were  taken  to  be  the  results  at  £  =  0.  At  f  =  2A£,  a  linear  extrapolation 
of  the  results  at  £  =  0  and  i  =  A{  was  used.  For  £  >  3A£  a  three  point  inverse  polynomial 
extrapolation  was  used  [see,  for  example,  Traub  ( 1 964) J.  This  approach  was  found  to  improve 
the  convergence  of  solutions  up  to  40%  faster  than  simply  using  the  previous  node's  values  as 
the  first  guess. 

To  improve  convergence  of  solutions  it  was  found  that  the  value  for  >\.r  needed  to  be¬ 
taken  larger  as  the  value  of  the  curvature  parameter  £  increase.  The  value  of  »/„,w as  initially  set 
to  10  for  Pr  >  0.7  and  to  15  for  Pr  =  0.1.  It  was  increased  by  5  when  the  value  of  P  or  </>  at 
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r\  =  0.98^  was  greater  than  0.001.  Since  P(^,oo)  =  <}>(£,  oo)  =  0,  the  P  and  <j>  values  greater 
than  0.001  at  y\  =  0.98*1^  indicate  that  the  corresponding  boundary  layer  thickness  extends  past 
>/,*,.  To  make  this  adjustment,  all  values  past  rj^were  set  to  zero  except  for  f  which  was  taken  to 
be  the  value  at  the  old  ^extended  out  to  the  new  >joo.  1'his  adjustment  introduced  less  than  0.5 
%  difference  in  the  values  of  the  local  Nusselt  number  when  compared  to  a  constant  run  at 
»7«,  =  30.  This  error  was  due  primarily  to  t initially  being  less  than  30  rather  than  to  the 
adjustment  itself. 

The  solution  method  is  an  iterative  scheme  and  a  solution  was  considered  to  be 
convergent  when  the  calculated  values  for  f,  P,  and  <f>  differed  from  the  last  guess  of  the 
respective  values  by  less  than  10  4  at  all  nodes  (i.e.,  at  all  r;  values  for  a  given  £).  When  these 
criteria  failed,  new  guesses  for  f,  P,  and  <f>  were  found  using  a  weighted  average  of  the  last  guess 
and  the  resulting  calculation.  That  is. 


fnew  =  (o  f  +  ( 1  -  «))fold 


f'«KW  =  ®P+(l-«)f' 


old 


4»new  =  «><£  +  (1  -  w)^>0|d 


(B  -  32) 


(B  -  33) 


(B  -  34) 


where  w  is  a  relaxation  factor.  Generally  w  =  1  resulted  in  quick  convergence.  I  lowever,  if 
was  to°  small  or  for  some  small  values  of  Pr  and  n  it  was  sometimes  necessary  to  use 
under- relaxation  and  set  w  =  0.5  to  facilitate  convergence. 

It  was  found  that  as  £  was  increased,  errors  resulted  due  to  increasing  boundary  layer 
thicknesses.  By  using  a  step  size  of  Arj  =  0.01  errors  were  reduced  for  calculations  at  t]  =  Oat 
high  values  of  the  curvature  parameter  £.  It  was  also  found  that  the  solution  was  not  sensitive 
to  the  step  size  in  £  and  Af  =  0.1  was  used. 
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APPENDIX  C 


TRANSFORMATION  OF  THE  GOVERNING  SYSTEM  OF  EQUATIONS  FOR  Till-; 
SPECIAL  CASE  OF  NATURAL  CONVECTION  IN  CHAPTER  3 


In  Chapter  3  the  solution  for  the  case  of  natural  convection  along  a  slender  vertical 
cylinder  with  variable  surface  temperature  is  presented  as  a  special  case.  Consider  a 
semi-infinite,  vertical  cylinder  with  radius  r0  that  is  aligned  in  a  quiescent  ambient  fluid  at 
temperature  T^.  The  axial  coordinate  x  is  measured  upward  when  Tw  >  and  downward 
when  Tw  <  TM.  Idle  radial  coordinate  r  is  measured  from  the  axis  of  the  cylinder.  The  surface 
of  the  cylinder  is  subjected  to  an  arbitrary  variation  in  temperature  T„(x),  and  the  gravitational 
acceleration  g  is  acting  downward.  Fluid  properties  are  assumed  to  be  constant  except  for 
variations  in  density  which  induce  the  buoyancy  force.  By  employing  the  laminar  boundary 
layer  assumptions  and  making  use  of  the  Boussinesq  approximation  the  governing  conservation 
equations  can  be  written  as: 

Continuity: 

Sr(n‘)+ 0  (C-‘> 


Momentum: 


—  d_/  rj3u\ 
r  dr  V  dr  / 


+  X/?(T-T00) 


(C~  2) 


Energy: 


M. 

dx 


4-  v 


dT 

dr 


(('-  3) 


In  these  equations  u  and  v  are  the  velocity  components  in  the  x  and  r  directions 
respectively;  T  is  the  fluid  temperature;  and  v,  and  a  are,  respectively,  the  kinematic  viscosity, 
the  volumetric  coefficient  of  thermal  expansion  and  the  thermal  diffusivity  of  the  fluid. 
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The  boundary  conditions  are 


u  =  v  =  0  ,  T  =  Tw(x)  at  r  =  rQ 


u  -*  0  ,  T  -» As  r  -*  oo 


u  =  0  ,  T  =  At  x  =  0,  r  >  r0 


(C-  4) 


(C-5) 


(C  6) 


In  writing  equation  (C-6)  it  is  assumed  that  the  flow  and  thermal  boundary  layer 
thicknesses  are  zero  at  the  leading  edge  of  the  cylinder  surface. 

The  conservation  equations  and  the  boundary  conditions  arc  then  transformed  into  a 
dimensionless  form  by  introducing  the  following  dimensionless  variables: 


1  =  -r°]  (Grx/4)1/4.  =  TL(Crx/4)-1  /4 


(C  -  7) 


f(T  ,)  =  ,i(x,r)/[4vr0(Grx/4)1/4J.  0(A,  n)  =  (T  -  T^)/[Tlv(x)  -  T^j  (C  -  8) 


Grx=J?/*[Tw(x)-Tw]  x3/*2 


(C-9) 


where  rj  is  the  pscudo-similartiy  variable,  f(A,  rj)  is  the  reduced  stream  function,  0(A,  >/)  is  the 
dimensionless  temperature,  A  is  the  curvature  parameter,  and  i^fx,  r)  is  the  stream  function  that 
satisfies  the  continuity  equation,  with  u  =  (cty/dr)/r  and  v  =  -  (di/f/dx)/r.  The  transformation 


u  =  TT(GrJt/4)l/2f'(4,  r\) 


(C  -  1«) 


'  -  -  f  f - ^)«  1 -  - 


(C  11) 


m 


L**V| 


+  _L  -.,^T  _  I'  )!/2~l  4-  xI|,2(T  -T  )i/2  f"—  +  — ^  — 

+  .X  («W  loo)  J  +  X  l‘w  *oo)  dx  +  ^  ,lx 


5© 

w 

Bt 
mm 


—  =  4v(Grx/4)3^4x~2-~-f"(,i,  £) 

V  lO 


8v(GrJC/4)3/V2-^f"  +  4v(Grx/4)x  3-^-P" 

2 


•'T  =(Tw(x)-T 


■>> 


(  ‘  w  -  *  J 


m  which  the  primes  denote  partial  differentiation  with  respect  to  rj. 


'VS.'T. 


Substitution  of  equations  (C-IO)  through  (C-18)  into  equations  (C-l)  through  ((2-6) 
results  in 

(1  4- "  +  Af '  +  (y  +  3)ff"  -  2(y  +  1)P2  +  0  =  A(y  -  l)(p'-|t-  -  (C  -  19) 

\  OA  OA  / 


(1  +  r,A)0"  +  ).0'  +  Pr(y  +  3)10'  -  Pr(4y)rfl  =  PrA(y  -  1)(»'-^  p|y) 


((•’-  20) 


iia  (»*  4*.  (4|  t  <  I.1 


1 


f(A,  o)  =  r(/L,  o)  =  o,  0(2,0)  =i 

P(2,oo)  =  0,  0(2,oo)  =  <) 


(C-  21) 


where 


(C-  22) 


The  system  of  equations  (C-19)  through  (C-22)  represents  the  general  form  of  the 
transformed  boundary  layer  equations  for  variable  wall  temperature  Tw(x)  along  vertical 
cylinders  in  natural  convection.  For  the  case  of  power  law  variation  in  surface  temperature, 


Tw(x)  =  T„  +  ax"  ,  one  obtains  from  equation  (C-22)  that 


(C  -  23) 


For  long  slender  cylinders,  the  curvature  parameter  2  =  2-^-(Grt/4)  1/4  can  be  large.  To 

*<) 

lower  the  maximum  value  of  calculations  in  the  x  or  2  coordinate,  one  introduces  a  new  £(x) 


variable  defined  by 


2  =  C£2 


(C  -  24) 


with  C  =  2«  4,/4. 


Substituting  equations  (C-23)  and  (C-24)  into  equations  (C-19)  through  (C-22)  results  in: 
Momentum: 

(1  +  a,»})f"  +  a,P'  +  a2fT  4-  a3f2  +  a40  =  a5(V-|^  -  f'-|^  (C  -  25) 


Energy: 


(1  +  a \r,)0"  +  a ,0'  +  Pr  a2f0'  +  Pr  a hf0  =  Pr  a 5^'-^  - 


Boundary  Conditions  : 


(C  -  26) 


Pr 

m 


& 


.Vc'' 

j>Ts 


'M 


v  vv  v'.-wa:.  ,  , , 


T.V7-J,  - 


where 


u{,o)  =  r«,o)  =  o.  o({,o)=i 

P«,oo)  =  0({,oo)  =  0 


aj  =  C<r  ,  a2  =  n  +  3  ,  a3  =  -  2(n  +  I) 
a4  =  1  ,  a5  =  £(n  -  l)/2  ,  a6  =  -4n 


(C-  27) 


(C  -  28) 


The  physical  quantities  of  interest  arc  the  local  Nussclt  number  N'u,  the  average  Nusselt 
number  NuL,  the  local  wall  shear  stress  tw,  the  axial  velocity  distribution  u,  and  the  temperature 
profile  0(£,  tj).  From  the  definition  of  the  local  Nussclt  number 


v  _  _Ax_  _  *-lw  x_ 

‘  x  k  \\-Tx  k 


(C  -  29) 


►  .f 


m 


along  with  qw  =  -  k(d'['ldr)„  ,  one  finds 


Nux(Gr</4)-,/4  =  -0'(£tO) 


The  average  Nusselt  number  is  obtained  from  the  expression 


Nu‘--TT‘tI 


Substituting  h  into  equation  (C-31)  and  carrying  out  the  integration,  one  finds: 


Nu,  (Gtl/4)_1/4  =  -  ir(3+n)/4  fLx(3+n)/V({,  ())x~  'dx 

J» 


From  equations  (C-7)  and  (C-24)  one  can  write 


(x/r0)'  n1i/8 


(C  —  30) 


(C-31) 


(C  *  32) 


(C-33) 


>r 


t 


m 


.-,v  r. 


• 

v  w 

VW 

v:^v 

»> 


.  ,811/(1  -n) 


r()[Gro^] 


(C  34) 


where  Gr0  =  Grx  at  x  =  r0;  that  is,  Gru  =  g(}(  a  r0n )  rj/v2.  Substituting  equations  (C-34)  and 
(C-35)  into  equation  (('-32),  one  gets 


NuL(GrL/4)-1/4  =  -  -1— it  0)d{ 

t  n  J0 


where 


iu=i  at  x=L, 

M  =  (5  +  3n)/(l -n),  N  =  (6+ 2n)/(n  -  1) 


As  £  -*  0,  equation  (C-36)  approaches 

N^L(GrL/4)-I/4^0  =  -  y^-0'(0,  0) 


Next,  the  local  wall  shear  stress  is  found  from 


Substituting  equation  (C-13)  into  equation  (C-39)  one  obtains 

Tw  =  ^«rx/4)3,4r({t0) 

X 

From  equation  (  C-10  ),  the  axial  velocity  distribution  can  be  written  as 

-Hr  =  4(Grx/4)!/2r(^,  r,) 

and  the  temperature  profile  is  given  by  0(£,  ij)  =  (T  -  TJift !„  -  T,J. 


(C-  36) 


(C  -  37) 


(C  -  3S) 


(C-39) 


(C  -  40) 


(C  41) 


APPENDIX  D 


TRANSFORMATION  OF  THE  GOVERNING  SYSTEM  OF  EQUATIONS  FOR  MIXED 

CONVECTION  IN  CHAPTER  3 

In  Chapter  3  the  solution  for  the  case  of  mixed  convection  along  a  slender  vertical  cylinder 
with  variable  surface  temperature  is  presented.  Consider  a  semi-infinite,  vertical  cylinder  with 
radius  r0  that  is  aligned  parallel  to  a  uniform,  laminar  free  stream  with  velocity  u  „  and 
temperature  T„.  The  axial  coordinate  x  is  measured  in  the  direction  of  the  forced  flow  and  the 
radial  coordinate  r  is  measured  from  the  axis  of  the  cylinder.  The  surface  of  the  cylinder  is 
subjected  to  an  arbitrary  variation  in  temperature  Tw(x),  and  the  gravitational  acceleration  g  is 
acting  downward.  Fluid  properties  are  assumed  to  be  constant  except  for  variations  in  density 
which  induce  the  buoyancy  force.  By  employing  the  laminar  boundary  layer  assumptions  and 
making  use  of  the  Boussinesq  approximation  the  governing  conservation  equations  can  be 
written  as: 

Continuity: 

t — (ru)  +  (rv)  =  0  (/>-!) 

dx  dr 


Momentum: 


V_ d_(  du\ 

r  dr  \  dr  / 


±gpc I' -To.) 


(/)  -  2) 


Energy: 


,*L 

dx 


+  vf- 

dr 


(D  -  3) 


The  positive  sign  in  equation  ( ID-2)  applies  to  upward  forced  How  and  the  negative  sign  to 
downward  forced  flow.  In  these  equations  u  and  v  are  the  velocity  components  in  the  x  and  r 
directions  respectively;  T  is  the  fluid  temperature;  and  v,  /?  and  a  are,  respectively,  the  kinematic 
viscosity,  the  volumetric  coefficient  of  thermal  expansion  and  the  thermal  dilfusivity  of  the  fluid. 


puuuuuum^  .luuwJWJiRiuuw  vw  w  ^  m  n."  ll"  towtto  TTOrrr^?  a«' 


'jr'jTUT'J  w'Lr'iyi'rJ  Turitrvm  *v  y-i ' 
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The  boundary  conditions  are 

u  =  v  =  0  ,  T  =  Tw(x)  at  r  =  rQ 

u  — »  0  ,  T-  T^  As  r  -*  oo 
u  =  0  ,  T  =  T^  At  x  =  0,  r  >  r0 


(/)  -  4) 


C>  -  5) 


(!)  -  6) 


In  writing  equation  (D-6)  it  is  assumed  that  the  flow  and  thermal  boundary  layer  thicknesses  are 
zero  at  the  leading  edge  of  the  cylinder  surface. 

The  conservation  equations  and  the  boundary  conditions  are  then  transformed  into  ■> 
dimensionless  form  by  introducing  the  following  dimensionless  variables: 


i  2  2\ 

lr  -  T°’  /D,J/ 2  ,  n  >/4\  ,  x 

,,  =  -^<Rc‘  +  t,r>  >• 


(/)  -  7) 


f(z,  »,)  =  ^(x,r)/[vr0(Rei/2  +  GrJ/4)],  0(7.,  n)  =  (T  -  T^IO'Jk)  -  T^]  (I)  -  8) 

“x  =  -^-.  x  =  (i  +  tti/4r'  (;)-g) 

Re2 

Grx  =  gft T»W  -  Tw]  x3/v2 ,  Rex  =  uMx/.  (/->  -  10) 

where  is  the  pseudo-similarity  variable,  z  is  the  dimensionless  axial  coordinate,  f(/., »/)  is  the 
reduced  stream  function,  0(z,r\)  is  the  dimensionless  temperature,  i^(x,  r)  is  the  stream  function 
that  satisfies  the  continuity  equation,  with  u  =  (cty/drj/r  and  v  =  —  (di/</dx)/r,  fl,  is  the 
buoyancy  parameter  which  varies  from  zero  for  pure  forced  convection  to  infinity  for  pure  free 
convection,  and  y  is  the  mixed  convection  parameter  which  varies  from  zero  for  pure  free 
convection  to  one  for  pure  forced  convection. 

The  transformation  yields 

u  =  uooX~  2f,(z.  rc) 


(/>  ID 


F’jr’jr’jrrjt-’jT’jrrjt  -v» 


(1  +  >/A)P"  +  AP'  +  -L(2  +  n  -  X)(y  +  l))fP' 

4 

-  y(i  -  x)(y  +  1)P2  ±  (1  -  xfo  =  -  - 1’^'  ) 


(1  +  Y]\)0"  -f  A  O'  +  ±L(2  +  (1  -  y)(y  +  1  ))fO'  -  PryPfl 


(/)-  20) 


=  —  I>r  y.(o'-^-  — 

\  ()/.  i)y.  / 


f(z,  0)  =  P(z,  0)  =  0,  0(/„  0)  = 

P(z,  oo)  =  x^<  0('/.,  oo)  =  () 


(/>  —  21) 


(/)  -  22) 


where 


A  =  2-f  (Rei/2  +  GrJ'4)- 

4  r\ 


(/>  -  22) 


is  the  surface  curvature  parameter  and 


V  Tw(x)X-  dx[lw(X)  {~] 


(/)  -  24) 


The  plus  and  minus  signs  in  front  of  the  term  (1  •  xY^  *n  cqaa:<k,n  (12  20)  now  represent 
buoyancy  assisting  flow  and  buoyancy  opposmg  flow,  respectively,  l  or  pure  forced  convection 
X  =  1  and  for  pure  free  convection  x  —  0- 

The  system  of  equations  (D-20)  through  (D-24)  represents  the  general  form  of  the 
transformed  boundary  layer  equations  for  variable  wall  temperature  Tw(x)  along  vertical 
cylinders  in  mixed  convection.  l  or  the  case  of  power  law  variation  in  the  surface  temperature, 
Tw(x)  =  4-  uxn,  one  has  from  equation  (D-24)  that 


liquations  (D-20)  through  (D-22)  contain  three  x-dependent  parameters,  /(x), 
A(x),  and  y(x).  These  parameters  can  be  related  to  a  single  x-dependent  parameter  s(x)  defined 


v.Sy> 

*  %'  v*J 

V'/.V.*, 


w 


K<< 

■VV’ 

,-W 

>  Va 


Thus,  one  has 


A  =  2t\ 


x  =  [i+^l+nr' 


where 


^0  = 


(ir„ 


Re 


l-n 


L  ■"•o  J 


1/4 


(/>-  29) 


with  Gre  =  Gr,  and  Re0  =  Re,  for  x  =  r0  Also,  the  right-hand  sides  of  equations  (D-20)  and 
(D-21)  become,  respectively 

4  \  di  di  ) 


Rewriting  equations  (D-20)  and  (D-21)  one  has: 

Momentum: 

(1  +a^)P"  +  a1f'  +  a2fr'+a3f2-fa4(?  =  a5^r'-|-_p|.^  (/)  22) 

Energy 

(H-a^)0"+a)f?'  +  Pra2W'+!,ra6P(7  =  Fra5^'-^  P-|^)  (/>  W 


Boundary  Conditions  : 


f(4,  0)  =  P(4.  0)  =  0,  0(4,0)=! 

f(4,oo)  =  a7  ,  0(4,  oo)  =  0 
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t 


l 

■i 

V 

s 

A 


(/)  34) 


where 


a,=A,  a2  =  i-(2  +  (l -jr)(n+  1))  ,  a3  =  ±(x  -  l)(n  +  1) 
a4  =  ±(l-x)4-  as  =  —  4/4  ,  ab  =  ~n,  a7  =  y2 


(!)  -  VS> 


With  A  and  y  related  to  4,  the  functions  f  and  0  in  equation  (D-32)  to  (D-34)  arc- 
functions  of  (4,  >;)  and  depend  on  three  constant  parameters  n,  Pr,  and  These  equations  arc- 
now  in  the  form  that  can  be  solved  by  the  method  proposed  by  Ire  et  al.  [6],  See  Appendix  0 
for  the  solution  method  as  applied  to  equations  (I)-32)  through  (D-34). 

The  physical  quantities  of  interest  include  the  local  and  average  Nussclt  numbers,  the  local 
and  average  friction  factors,  the  axial  velocity  distribution,  and  the  temperature  profile.  I  lie- 
local  Nusselt  number  is  defined  by 

Nux  =  ^  (/>  -  36) 

where  h  =  qw/(Tw  T^)  and  q„  =  -  k(<7\'lih)t  t  .  Thus,  it  can  be  shown  that 

Nux/(RcJ/2  +  C.r'J*)  =  -  d'(4,  0)  (/>  -  37) 

The  average  Nussclt  number  is  obtained  from  the  expression 

Nu,  =  %  (/>  -  38) 

k 

where 

h  =  - --  (" '  /id X  (/>-  34) 

I  J0 


This  results  in 


Lt'LlVl'*.Vl 


^L-f, 

Jn 


Re^V#'^,  0)x_ldx 


NV(Rc|«  +  an  -  -  n.!-'^^-'. 


where  Xl  =  X  at  x  =  L.  F;rom  equation  (D-26) 


x  =  r„Re0£4 


x_ldx  =  4<r'd£ 


Substituting  equations  (D-42)  and  (D-43)  into  equation  (D-41),  one  arrives  at 


where  <*;L  =  £  at  x  =  L.  As  <*  -+  0,  equation  (D-44)  approaches 


NuL/(Re[/2  +  GrJ/4)|_0  =  -20'(O, 0) 


The  friction  factor  is  obtained  from  the  definition 


2tw  _  2  du  I 

P*lo  r=r° 


Substituting  equation  (D-14)  into  equation  (D-46)  results  in 


Cf ReJ/2  =  2x“3f ’((,  0) 


The  expression  for  the  average  friction  factor  is  derived  from 


h.  2 
P»oo 


_rlTwdx=^rl^l  ax 

I,  Jo  uLl-Jo  dT  'r=rn 


(/)  -  40) 


(D-  41) 


(/>-  42) 


(D  -  43) 


Nu  '(Re[/2  +  Gr^4)  =  -4xL^2f ' '  (/>  -  44) 


(/)  -  45) 


(/)  -  46) 


(/)  -  47) 


(/)  -  48) 
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Substituting  equation  (D-14)  into  equation  (D-48),  one  obtains 


Cf  Rc[/2  =  2L_I/2  l'  xl2x-3["(t,  0)x_ldx 

J(\ 


(/)-  49) 


Next,  substituting  equations  (D-42)  and  (D-43)  into  equation  (D-49),  one  obtains 


CfLRel/2  =  8r2J’V3n£.0)£d£ 


(D-  50) 


As  {  -*  0,  equation  (D-50)  approaches 


Cf1Rci./2^o*4HO,0) 


From  equation  (D- 11),  the  axial  velocity  distribution  can  be  written  as 

U  fU  r,) 


and  the  temperature  profile  is  given  by  0(£,  tj)  =  (T  -  T00)/(TW  - 


(D  -  51) 


(l>  -  52) 


ft 

£■ 

•&T 


S’.O 


I 


WSw 


9t 


5V 

M 

,*vv- 


'  v> ' 
<nV- 


% 

M 


APPENDIX  E 


SOLUTIONS  ()!•  THU  SYSTEM  OI  EQUATIONS  (C-25)-(C-27)  AND  (D-32)-(D-34) 


Equations  (C-25)  through  (C-27)  are  identical  in  form  to  equations  (D-32)  through 
(D-34),  with  the  coefficients  a,  through  a6  given  by  (C- 28)  for  pure  free  convection  and  by 
equation  (D-35)  for  mixed  convection.  The  solution  method  of  Ire  et  al.  [6]  was  employed  to 
solve  the  two  systems  of  equations. 
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: "\>vyv 
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<•  V ' 
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The  first  step  is  to  convert  the  terms  involving  in  the  following  manner. 


~  =  (H  —  II0)  —  q(-^ii-)0  =  — (H  —  H  ) 

d{  A<*  0  d{  °  A£  v  0/ 


U:--  1) 


<V${ 

"W  N.  * 


where  II  is  a  function  of  (£,  >j),  II0  =  II0  +  (q/p)A£(dII/d£)0,  and  the  subscript  "o"  denotes 
quantities  at  £  -  A£.  The  values  for  p  and  q  arc  p~  1  and  q  =  0  al  £  =  0  and  £  =  A£, 


and  p  =  2  and  q  =  1  thereafter  for  <*  >  2A£.  Thus  we  have 


A{  o) 


(T- 2) 
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rvyv-rv 

)•  V*  v*  V 
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—  =  — (0  -  0O) 


(T  -  3) 


(T  -  4) 
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0  =  0  +  —  M(—) 

0  0  +  p  Q{  d£  0 


(/' -  7) 


The  f  and  0  equations  are  then  linearized  by  using  the  simple  technique  for  the  product  of 
two  arbitrary  functions  F  and  G  ,  defined  as 


FG  =  FG  +  GF  -  FG 


(/-'  -  X) 


where  F  and  G  indicate  the  values  to  be  guessed.  Applying  equation  (F-8)  to  the  terms 


ff”  ,  P2 ,  f O' ,  and  P 0  ,  one  obtains 


flf"  =  ff"  +  f"f-  f"f 


P2  =  2PP  -  P2 


W  =  {  O'  +  O' f-  f O' 


ro  =  ro  +  ro-ro 


(A  -  9) 


(/■:-  id) 


(A  -  ii) 


(A  —  12) 


Substituting  equations  (F-2)  through  (F-12)  into  equations  (C-25)  and  (C-26)  (or  I >-32 
and  D-33)  results  in: 

Momentum  equation: 

(1  +  *,»f)P"  +  Cai  +  (a2  -  as  P/Ai)  f  +  a5  p/A£  fQ]  f" 

+  [2(a3  +  a5  p/Af)  P  -  a5  p/A{  P0]  P  +  [(a2  -  a5  p/A£)f']  f  +  a40  (/:  -  13) 

=  C(a2  -  as  P/A£>fP'  +  (a3  +■  a5  p/A£)f ,2] 


Energy  equation: 


(1  +  a,»j)<£"  +  [a,  +  l’r(a2  -  asp/A£)f  +  Pr(a5p/A0fo]fP 
+  [Pr(a6  +  a5p/Af  jf']0  +  [IT{(a6  +  a5p/A£)0  -  a5p/A{0o}]P 
+  [Pr(a2  -  asP/A$)0']f  =  [Pr{(a2  -  asp/A{jfO'  +  (a6  +  asP/A£)f'0}] 


(A  -  14) 
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Equations  (E-13)  and  (E-14)  can  be  rewritten  in  the  following  form. 


A0P"  +  (a,  +  A,)f"  +  A2P  +  A3f  +  A40  =  A5  (A  -  15) 


iyr  +  (a,  +  n,)0'  +  b2o  +  n/  +  B4f  ==  b5  a:  i^ 

where 

A  q  —  l  4-  a  1 1/ ,  A|~ a+f  +  d*f0 

A2  =  2b+f'  -  d*Pu  ,  Aj  =  a*P'  {E  -  17) 

A4  =  a4  ,  A5  =  a*fP'  +  b*f'2 

B0  =  A0,  B,  =  Pr  A,  ,  B2  =  Pr  c*f' 

B3  =  Pr  (c *0  -  d *0o) ,  B4  =  Pr  a *0'  (E  -  IK) 

B5  =  Pr(a*f0'  +  e*f'0) 


a*  —  a2  —  d* ,  b*  =  a3  4-  d* 
d*  =  a5p/A^,  c*  =  a6  +  d* 


(/-  -  19) 


The  coefficients  a,  through  a6  and  the  functions  f0,  f„.,  and  0U  are  as  previously  defined. 

The  f  and  0  equations  have  now  been  reduced  to  a  system  of  quasi-linear  ordinary 
differential  equations.  By  defining  g  =  P  and  applying  the  weighting  factors  (see  I,ec  et  al.  [6J), 
one  can  write  equations  (E-15)  and  (E-16)  in  the  finite-difference  form  as: 

f1-fi._,-(A»,/2)(gi  +  gI._1)  =  ()  (A  -  20) 

Ao/(Afj)2(a-lg,-l  +  aog<  +  “-t-lBi-l-l)  +  ^2g,  +  A3f(- +  A40j  =  A5  (A  -21) 

BO/(A»7)2(/?_|0;_|  +  Mi  +  /?+i^+1)  +  B20,  +■  B3&  +  B4f(  =,  B5  (/•;  -  22) 


where  Aq  is  the  step  size  in  the  tj  direction  and  the  subscript  "i"  refers  to  values  at  the  nodal 
point 


In  equations  (E-21)  and  (E-22)  the  weighting  factors  are  defined  as: 

o+1  =  W,<ZI+|/2),  a0  =  -a_,-a+1 

P-i  -  WK  -  /  V 1/2)  ■  P+\  =  wK*Vi/2>  ■  ft.  =  -  P-i  -  P+\ 

W({z)  =  z/(l  -cxp(-z)) 
z  =  A>f(a,  +  A,)/A0  ,  z*  =  A>/(a,  +  B,)/B0 


(/-  -  23) 


Equations  (E-20)  through  (E-23)  are  applied  to  the  interior  points.  On  the  boundaries, 
equation  (C-27)  or  (D-34)  applies  such  that 


f«=R.  =0,  0,=  1 


fn-fn-l  --T-(gn  +  8n-l)  =  °n  =  l).  gn  = 


(A- -24) 


where  the  subscript  n  is  the  number  of  nodes  in  the  >7  direction,  and  a-  =  0  for  natural 
convection  and  a,  =  x2  for  mixed  convection. 

Equations  (E-20)  through  (E-24)  constitute  a  system  of  algebraic  equations  that  can  be 
written  in  the  matrix  form 


[A][X]  =  [B] 


(I-  25) 


where  [A]  is  a  band  matrix  of  order  3n  and  bandwidth  seven.  The  array  [X]  which  contains 
the  solution  in  the  form  (f,,  g, ,  0,,  f2,  g2,  02,  ....  fn,  g„,  0„)T  is  a  column  matrix  of  order  3n. 
Idle  matrix  [B]  is  a  column  matrix  of  order  3n  which  contains  the  right  hand  sides  of  equations 
(E-20)  through  (E-22)  and  (E-24).  The  matrix  [A]  is  approximately  diagonally  dominant  and 
equation  (E-25)  can  be  solved  by  the  Gaussian  elimination  technique  with  high  accuracy. 

To  obtain  values  for  f"  and  O’,  the  results  for  P  and  0  along  with  the  boundary  conditions 


P”(£,  0)  =  —  a,P'(£,  0)  —  a4 


0"«,O)  =  -ai0UO) 


P"(£,  00)  =  <?"(£,  00)  =  0 


(/•:  -  26) 


(/•;  27) 


(/■:  -  2K) 


are  used  in  a  cubic  spline  interpolation  routine  (see,  for  example,  Burden  and  1'aires  [7]  ).  As 
the  solution  converges,  the  boundarv  conditions  become  exart  and  the  values  for  P'  and  O'  „.m 
be  obtained  with  high  accuracy. 


The  present  method  of  solution  employs  a  quasi-linearization  of  the  original  nonlinear 
system  of  equations  and  requires  initial  guesses  for  f,  P,  P',  0,  and  O'.  The  flat  plate  solution 
(i.e.,  £  =  0)  for  the  uniform  wall  temperature  case  (L’VVT)  for  Pr  =  0.7  was  obtained  and  these 
results  were  used  as  the  initial  guesses  for  all  other  combinations  of  Pr,  n,  and  fi0  at  f  -  0.  I  or 
>  A£,  good  convergence  was  obtained  by  simply  letting  the  solution  at  the  previous  node  be 
the  guesses  for  the  next  node. 

The  solution  method  is  an  iterative  scheme.  A  solution  was  considered  to  be  convergent 
when  the  calculated  values  for  f,  P,  and  0  differed  from  the  last  guess  of  the  respective  values  by 
less  than  10  4  at  all  nodes.  When  these  criteria  failed  to  be  met,  new  guesses  for  f,  P,  and  0 
were  found  using  a  weighted  average  of  the  last  guess  and  the  resulting  calculation.  That  is 

fnew  =  o>  f  +  (1  -  w )f0|d  (/-  29) 

f'new  =  a,P+(l-c0)f'old  (!■:  -  30) 

0ncw  =  °>0  +  (\  -  (d)0o ld  (A.  -  ,11) 

Generally  <±>  =  1  resulted  in  quick  convergence.  However,  for  =  2  it  was  often  necessary  to 
use  under-relaxation  with  to  =  0.5  to  facilitate  convergence. 

It  was  found  that  numerical  results  depended  upon  the  choice  of  As  £  was  increased, 
it  was  necessary  to  increase  The  value  of  was  initially  set  to  1 5  and  was  increased  by  5 
when  the  value  of  |P  — a,|  or  0  was  greater  than  0.001  at  rj  =  0.98^.  Since 
P(£,  oo)  —  a7  =  6(£,  oo)  =  0,  the  |  P  —  a7 1  and  0  values  greater  than  0.001  at  0.98  indicate  that 
the  corresponding  boundary  layer  thickness  extends  past  rj^.  To  make  this  adjustment  for 
natural  convection,  all  values  of  f  and  0  and  their  derivatives  past  ^  were  set  to  zero  except  for 
f  which  was  taken  to  be  the  value  at  the  old  extended  out  to  the  new  »}„. 


mpuww.irwvi  F.-ivisT; 


The  adjustment  for  mixed  convection  is  slightly  more  complicated.  Again  all  values  were 
set  to  7er o  except  the  values  for  f,  P,  (df[df)0,  f  and  i*  .  They  were  adjusted  as 

follows. 


f,.  =  f,.1+A^V 


fwv. 


(I)--'7'- 


(jy)°i  =  -2xoii°(n+  IKo 

foi  =  foi  +  A£  ( ■ oi 


fo  i  +  -s-A{ 


(f> 


(i;  32) 


(/:’  —  33) 


(/-  -  34) 


(/■:  -  35) 


(/'  -  36) 


(/•:  -  37) 


This  adjustment  introduced  less  than  1.0  %  difference  in  the  local  Nusselt  number  for  l*r  =  0.7 
as  compared  to  a  constant  run  at  >7  =  30.  This  error  was  due  primarily  to  rj^  being  less  than  30 
initially  rather  than  to  the  adjustment  itself. 

It  was  found  that  as  £  was  increased,  errors  resulted  due  to  the  growing  boundary  layer 
thicknesses.  By  using  Atj  =  0.01  numerical  errors  were  reduced  for  calculations  at  high  values  of 
the  curvature  parameter  It  was  found  that  the  step  size  A<j;  did  not  affect  the  solution 
appreciably  and  A£  =  0.1  was  used. 
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